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Abstract 


This  report  is  directed  toward  the  design  of  a real- 
time estimation  algorithm,  a Kalman  filter,  that  estimates 
aircraft  position  and  velocity  using  multiple  DME  range 
measurements.  The  estimator  is  designed  and  tested  for 
feasibility  as  a reference  system  for  examining  Inertial 
Navigational  System  (INS)  low  frequency  errors.  Both  a 
9-state  estimator  including  jerk  states  and  a 7-state  esti- 
mator without  the  jerk  states  are  designed. 

With  the  tuning  parameters  used  in  the  estimator  tests, 
the  7-state  estimator  provides  better  performance  than  the 
9-state  estimator.  An  approximate  analysis  of  the  7-state 
estimator  performance  (by  comparison  to  FASTMAP,  a currently 
used  and  accepted  filter,  and  CIRIS,  the  Completely  Inte- 
grated Reference  Instrumentation  System) , reveals  that  esti- 
mator errors  in  the  high  frequency  range  are  greater  than 
those  of  an  INS,  but  errors  in  the  DME-based  estimator  are 
consistent  in  strength  and  do  not  exhibit  an  unbounded 
growth  as  typical  of  INS  errors.  For  the  estimator  in  this 
study,  the  approximate  values  that  encompass  50  percent  of 
all  the  errors  (as  compared  to  CIRIS)  for  latitude,  longitude, 
latitude  velocity,  and  longitude  velocity  were 
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Latitude  position 
Longitude  position 
Latitude  velocity 
Longitude  velocity 


+ 83  feet 
+183  feet 
+ 8,4  feet/sec 

+ 7.5  feet/sec 


Improving  estimator  performance  is  suggested  by  proper 
tuning  and  by  using  an  adaptive  approach. 


A DESIGN  OF  A TRAJECTORY  ESTIMATOR 


USING  MULTIPLE  DME  RANGE 
MEASUREMENTS 

I.  Introduction 

An  important  component  of  the  modern  navigation  system 
acquisition  process  is  the  flight  test  program.  The  flight 
test  program  is  used  to  evaluate  and  verify  the  performance 
of  inertial  navigational  systems  (INS)  and  includes  the 
proper  selection  of  reference  navigational  systems.  Current 
flight  test  reference  systems  employ  externajyroeasur^Bments 
such  as  radar,  Distance  Measuring  Equipment /(DME)  .^corre- 
lation techniques,  and  onboard  reference  INSv-  These  meas- 
urements are  used  to  form  reference  trajectories  which  are 
compared  with  the  INS  trajectory  data  to  evaluate  the  INS. 

Depending  on  the  type  of  reference  system  used,  such  an 
analysis  can  either  be  a post-flight  evaluation  or  a real- 
time evaluation.  A real-time  evaluation  of  the  system  keeps 
the  pilot  (or  operator)  continuously  aware  of  the  system's 
performance.  In  this  way,  an  INS  malfunction  is  quickly 
detected,  and  perhaps  a costly  mission  is  aborted. 

Trajectory  errors  of  high  frequency  are  relatively 
uncommon  in  INS;  that  is,  INS  short-term  oscillations  are 
minimal.  In  contrast,  INS  trajectory  errors  of  low  frequency 
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are  of  substantial  importance.  (The  84-minute  Schuler  period 
is  always  evident  in  the  flight  tests.)  Although  inertial 
navigational  systems  may  be  highly  accurate  for  a short  time 
after  initialization  (alignment),  these  systems  are  hindered 
by  long-term  drift  errors  that  increase  with  time.  A refer- 
ence system  for  the  INS  should  be  more  accurate  than  the  INS 
in  the  low  frequency  error  domain  to  provide  an  adequate 
method  for  checking  INS  performance  for  typical  flight  time. 
Since  high  frequency  errors  in  the  INS  are  usually  insignif- 
icant, INS  superiority  to  the  reference  in  the  high  frequency 
range  is  not  intolerable. 

The  need  for  a low-cost,  real-time,  and  easily  deploy- 
able reference  system  (with  a relatively  small  low  frequency 
error)  has  led  to  investigating  the  use  of  existing  DME  sta- 
tions as  a source  of  reference  systems  (Refs  3:  4;  5).  DME 
stations  are  capable  of  giving  noise-corrupted  range  meas- 
urements to  aircraft  by  returning  a signal  received  from  the 
aircraft.  A continuous  input  of  local  station  ranges  can 
conceivably  be  employed  in  a minicomputer  algorithm  or  cen- 
tral processor  to  produce  aircraft  trajectory  estimates. 

This  multiple  DME  reference  system  would  have  small  low  fre- 
quency errors  since  DME  errors  are  rather  consistent  in  RMS 
magnitude  and  are  not  characterized  by  unbounded  error 
growth  as  typical  in  inertial  navigation  systems.  A refer- 
ence system  using  DME  stations  could  also  exploit  bearing 
data  for  a trajectory  determination.  However,  bearing  data, 
available  when  the  DME  stations  are  a part  of  VOR/DME  or 
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TACAN,  are  generally  much  noisier  than  range  data.  In  light 
of  the  noisy  nature  of  this  data,  a trade-off  exists  between 
keeping  the  reference  as  simple  a&  possible,  and  increasing 
complexity  by  the  addition  of  more  information.  In  consider- 
ation of  the  reference  system  in  this  study,  trajectory  infor- 
mation that  could  be  extracted  from  bearing  data  is  simply 
neglected  in  favor  of  avoiding  a more  complex  design.  Only 
DME  range  measurements  are  used. 

Available  local  station  ranges  include  measurements  from 
stations  within  150  nautical  miles  for  aircraft  attitudes 
above  18,000  feet  (Ref  5:150).  (Greater  ranges  can  be 
obtained  for  higher  altitudes . ) Every  point  in  the  United 
States,  except  for  portions  of  the  Northwest,  is  covered  by 
at  least  ten  DME  stations  within  a 150  nautical  mile  radius 
(Ref  5:150).  Because  over  750  DME  stations  already  exist 
in  the  continental  United  States  as  a part  of  VORTAC,  VOR/DME, 
and  TACAN  facilities  (Ref  l:Chap.  IX,  p.  8),  the  transition 
to  this  type  of  reference  system  should  be  reasonably  quick 
and  inexpensive. 

This  study  is  directed  toward  the  design  of  a real-time 
estimation  algorithm,  a Kalman  filter,  that  estimates  position 
and  velocity  of  an  aircraft  using  multiple  DME  range  meas- 
urements. In  addition  to  position  and  velocity,  the  algorithm 
also  estimates  a bias  associated  with  each  station  measurement. 
In  order  for  the  estimator  to  be  a feasible  reference  compared 
with  currently  used  references  (see  Chapter  II),  accuracy 
goals  are  set  to  encompass  half  of  the  longitude  and  half  of 


the  latitude  position  errors  within  +100  feet  and  half  of 
the  velocity  errors  within  +8  feet  per  second.  Such  accu- 
racy in  the  reference  should  verify  long-term  INS  errors; 
(typical  INS  can  have  one  nautical  mile/hour  drift  rates) . 
Although  the  Kalman  filter  is  designed  and  analyzed  in 
FORTRAN  on  the  60-bit  CDC  6600  at  Wright-Patterson  Air  Force 
Base,  an  extension  of  this  study  includes  the  writing  of  the 
algorithm  into  assembly  language.  Bearing  in  mind  the  word- 
length  problem  (numerical  precision  and  numerical  stability) 
associated  with  converting  from  60  bits  to  16  bits  (or  so), 
a minicomputer  or  a general  purpose  machine  already  onboard 
can  then  be  used  in-flight  to  process  the  incoming  range 
data  with  the  algorithm. 


II . Background 


History 

In  December  1976,  the  Central  Inertial  Guidance  Test 
Facility  (CIGTF)  conducted  flight  tests  of  the  FASTMAP  (Fast 
Multi-DME  Airborne  Position)  system.  The  FASTMAP  tests 
involved  the  use  of  multi-DME  range  measurements  to  compute 
aircraft  position  and  velocity.  FASTMAP  system  operations, 
once  initiated,  were  automatic;  the  DME  airborne  interrogator 
automatically  sequenced  frequencies  of  stations  in  the  vicin- 
ity of  the  flight.  A frequency  corresponds  to  the  identifi- 
cation number  of  a particular  station.  The  channel  number 
of  each  DME  station,  system  time,  signal  power  level,  atmo- 
spheric temperature  and  pressure,  and  each  noise- corrupted 
measurement  were  stored  in  a raw  data  package.  Trajectory 
data  from  CIGTF 's  Completely  Integrated  Reference  Instru- 
mentation System  (CIRIS)  were  also  stored  for  the  same  flights. 
The  position  and  velocity  computations  were  accomplished  post- 
flight and  compared  to  CIRIS.  Position  and  velocity  accuracy 
obtained  from  the  FASTMAP  system  were  109.4  feet  CEP  (Circular 
Error  Probable)  and  9.2  feet  per  second  CEP  respectively 
(Ref  3:21). 


Use  of  Multiple  DME  for  a Position  Determination 

An  actual  DME  measurement  involves  determining  the  time 
required  for  a radar  signal  to  travel  from  the  aircraft  to 
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the  DME  station  and  back  to  the  aircraft.  Knowledge  of  the 
signal  propagation  velocity  (C)  and  time  lapse  (At)  can  be 
used  for  a rough  distance  measurement.  The  common  formula 
D = CAt  is  employed  for  the  range  calculation.  Each  station 
has  a nominal  50  microsecond  delay  which  must  be  subtracted 
from  At  automatically  or  by  an  operator  in  the  aircraft 
(Ref  5:151). 

Although  many  range  measurements  should  be  available  to 
the  system,  the  type  of  equations  to  be  employed  in  the 
position  determination  can  be  illustrated  best  by  the  use  of 
three  DME  range  measurements  in  the  following  simplified 
example.  Without  the  presence  of  system  errors,  three  (or 
more)  DME  range  measurements  can  be  used  to  determine  an 
aircraft's  position  by  triangulation  methods.  This  require- 
ment of  triangulation  no  longer  applies  when  a dynamic  model 
is  introduced,  as  in  the  Kalman  filter.  In  other  words,  with 
a dynamic  model,  flight  trajectory  information  is  attainable 
from  DME  measurements  taken  one  measurement  at  a time . 
Nevertheless,  the  equation  required  for  each  of  the  three 
range  measurement  determinations  in  the  example  is  essen- 
tially the  same  equation  used  for  the  dynamic  model. 

In  Figure  1,  a diagram  for  the  acquisition  of  three  DME 
stations  is  shown,  €2,  and  are  the  actual  DME  meas- 

urements, r,j,  is  the  position  vector,  and  the  coordinate 
frame  xyz  is  arbitrary.  Three  equations  can  be  solved  for 
r»,  , r , and  r in  terms  of  e.  , e , and  e and  station 
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Figure  1.  Position  Via  Three  DME  Stations 


position  coordinates,  r.  , r_  , r_  : 

i ®i  ^i 


= (r_  -r.  + (r  -r.  + (r  -r 

j il  1.2  ^2  ^3  ^3 


The  subscripts  A,  B,  and  C represent  the  three  different  DME 
stations,  and  the  subscript  T represents  the  aircraft 
position;  i = 1,  2,  or  3 denotes  the  particular  component. 


In  practice  it  is  impossible  to  acquire  three  stations 
at  once.  Therefore,  the  algorithm  for  this  study  must  have 


the  capability  ,tQ . calculate  positions  from  a single  meas- 
urement a finite  time  apart  from  another  measurement.  Also, 
unlike  the  simplified  example,  the  algorithm  design  needs 
initial  positions  and  velocities  to  start  the  estimating 
procedure.  Though  the  approaches  are  quite  different,  a 
range  measurement  equation  similar  to  the  ones  used  in  the 
example  is  employed  in  the  algorithm  as  the  measurement 
equation. 

DME  Error  Model 

In  the  preceding  example,  measurement  errors  were 
assumed  absent.  However,  associated  with  each  DME  range 

K measurement  is  an  uncertainty  due  to  several  different  types 

V 

of  errors.  These  errors  can  be  separated  and  analyzed  to 
form  an  overall  DME  measurement  model.  The  following  DME 
error  model  is  credited  to  the  investigation  of  DME  errors 
by  R.  W.  Latham  and  R.  S.  Townes  in  1975  (Ref  6:332-342)  and 
others  (Refs  2;  7). 

The  error  in  each  DME  measurement  consists  of  errors  in 
the  airborne  equipment,  the  propagation  path,  and  the  ground 
station.  Latham  has  devised  an  error  model,  based  on  pre- 
vious models,  that  separates  each  of  these  error  types  into 
a bias  error  and  a wideband  noise  error.  Bias  errors  consist 
of  constant  errors  which  cause  the  range  measurement  to  be 
always  more  or  less  than  the  true  value.  In  contrast,  noise 
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errors  randomly  oscillate  about  the  bias  and  can  change  from 
measurement  to  measurement.  The  four  major  error  contri- 
butions have  been  found  by  past  experience  to  be:  bias 
errors  in  the  ground  station  (394  feet  RMS) , bias  errors  in 
the  airborne  equipment  (164  feet  RMS) , noise  in  the  airborne 
equipment  (50  feet  RMS) , and  noise  in  the  ground  station 
(26  feet  RMS)  (Ref  6:332). 

DME  ground  stations  are  intended  to  transmit  a radio 
signal  exactly  50  microseconds  after  receiving  a signal  from 
the  aircraft.  The  50  microsecond  delay  comes  from  natural 
delays  in  the  electronic  equipment,  a delay  line,  and  a 
finely  adjustable  electronic  delay.  Any  deviation  from  the 
50  microsecond  delay  will  cause  an  error  in  the  range  meas- 
urement. Inaccuracies  inherent  in  the  ground  equipment  are 
^ responsible  for  such  deviations. 

k’ 

Airborne  equipment  errors  are  caused  mainly  by  power 
level  uncertainties.  Latham  has  shown  that  as  the  power 
level  increases,  the  bias  errors  change.  Because  of  this 
functional  dependence  of  bias  on  signal  strength,  the  amount 
of  free  space  attenuation  also  affects  the  error. 

The  development  of  a DME  error  model  now  allows  the  use 
of  Kalman  filtering  to  calculate  the  best  estimate  of  the 
position  and  velocity  states.  Bias  errors  are  for  the  most 
part  removed  by  repeated  experiments  (see  Chapter  IV) , and 
the  mean  value  of  noise  errors  is  assumed  to  be  zero. 

Whereas  the  example  in  this  chapter  presents  a static 
and  noiseless  situation,  a Kalman  filter  must  take  into 
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account  the  dynamic  quality  of  an  aircraft  in  flight  and  an 
error  model  such  as  the  one  described  above.  The  Kalman 
filter  in  this  study  is  constructed  to  give  estimates  of 
the  trajectory  states.  These  estimates  are  essentially 
statistically  weighted  averages  of  the  solution  to  a set  of 
dynamic  equations  and  of  DME  range  measurement  information. 
The  algorithm,  or  Kalman  filter,  is  henceforth  referred  to 
as  an  "estimator." 


ff 


III . Estimator  Design 


Theory 

Now  that  a noise  model  and  measurement  equation' have  been 
developed  in  the  previous  chapter,  the  estimator  can  be  designed 
to  meet  the  prescribed  accuracy  goals.  Before  actually  building 
the  estimator,  the  dynamics  equations,  measurement  equation,  and 
statistical  characterization  of  noises  and  uncertainty  need  to 
be  specified.  The  estimator  uses  information  from  both  the 
dynamic  equations  and  DME  measurements  to  obtain  best  estimates 
of  the  aircraft  position  and  velocity  states. 

An  alternative  to  estimating  trajectory  states  is  to  estimate 
error  states,  such  as  INS  states  minus  estimator  states.  A 
reference  system  using  this  approach  can  conceivably  be  used 
for  evaluating  INS  output.  However,  the  reference  designed 
in  this  study  estimates  total  trajectory  states,  and  comparison 
to  INS  states  is  accomplished  outside  of  the  estimator  algorithm. 

For  the  most  part,  the  usual  methods  for  extended  Kalman 
filter  (EKF)  design  are  employed  to  keep  the  design  straight- 
forward. Nevertheless,  several  ad  hoc  procedures  are  necessary 
for  this  estimator  problem.  For  example,  the  measurement  com- 
putations require  an  evaluation  of  the  range  equation: 


e 


- ’'a  + (’^1 


- r, 


) + (r. 


- r, 


+ b.  + V 


(4) 


where  b^  is  the  DME  measurement  bias  associated  with  station  i 
and  V is  zero-mean  white  noise.  Use  of  Equation  (4)  requires 
keeping  track  of  which  station  is  being  acquired,  both  for  the 
appropriate  station  coordinates  and  the  bias  evaluation.  Since 
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the  bias  is  unique  for  each  station,  a procedure  is  necessary 
to  keep  the  bias  and  bias  variance  for  each  station  separate. 

Both  a 9 - state  estimator  and  a 7 - state  estimator  are 
designed.  The  synthesis  of  the  9 - state  estimator  includes  as 
states  the  first  three  derivatives  of  the  position  states: 
velocity,  acceleration,  and  rate  of  change  in  acceleration 
(jerk) . The  7 - state  estimator  includes  the  first  two 
derivatives  of  the  position  states:  velocity  and  acceleration. 
First,  the  estimator  that  includes  jerk  is  designed.  The  9 - 
state  estimator  is  then  easily  transformed  to  the  7 - state 
estimator  by  modelling  the  noise  as  entering  at  the  next 
lower  derivative  level.  (See  page  21).  Both  estimators  are 
tested  for  performance  and  the  final  choice  between  these  are 
made  in  Chapter  IV. 

In  every  estimator  problem,  a suitable  coordinate  frame 
and  appropriate  units  must  be  chosen.  Two  different  approaches 
are  investigated.  One  possible  coordinate  frame  for  the  estimator 
states  is  an  XYZ  eartesian  frame  centered  at  the  earth's  center 
with  the  Z - axis  through  the  north  pole,  the  Y - axis  through 
the  Greenwich  Meridian  and  equator,  and  the  X - axis  forming  a 
right  - handed  coordinate  system.  Position  inputs  would  be  in 
terms  of  latitude,  longitude,  and  height  from  the  local  geodetic 
frame  of  reference,  but  they  would  be  converted  to  the  XYZ  frame 
by  the  estimator  using  an  oblate  earth  model.  At  the  conclusion 
of  each  estimation  process,  the  updated  states  in  XYZ  coordinates 
would  be  converted  back  to  the  familiar  latitude,  longitude,  and 
height  for  output. 
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Another  possible  approach  is  to  skip  the  input/output 
conversion  required  above  and  perform  the  estimation  process 


directly  in  the  local  geodetic  frame  using  latitude  and  longitude 
in  the  estimator  equations.  For  this  study,  working  directly 
in  the  local  geodetic  frame  proves  to  be  most  useful  because  only 
one  conversion  from  the  local  geodetic  frame  to  the  XYZ  frame 
is  used.  The  es;  'cor  performs  entirely  in  the  local  geodetic 
frame,  but  one  conversion  is  needed  to  utilize  the  DME  information 
in  the  distance  relationship.  (See  page  25)  Table  I illustrates 
the  units  of  each  quantity  used  in  the  estimator.  The  last 
column  contains  actual  converted  output  units;  a blank  in. this 


Table  I.  Units  of  Estimator  States 


QUANTITY 


Angular  Position 
Angular  Velocity 
Angular  Acceleration 
Angular  Jerk 
Bias 

Position  Covariance 
Velocity  Covariance 
Acceleration  Covariance 
Jerk  Covariance 


ESTIMATOR 

UNITS 


radians 

radians/ second 
radians /second^ 
radians/ second^ 
radians- 

radians'^  2 

(radians/ second) 
(radians/second^) 2 
(radians/ second^) 2 


OUTPUT 

UNITS 


radians 
knots  or  feet 


feet 

feet*^  or  radians'^ 


column  signifies  there  is  no  output  for  the  listed  quantity. 

The  estimator  synthesis  necessitates  the  use  of  five 
basic  filter  equations: 

i (k+l)  " i ^tk) 


- (k+l) 


= <I>  + G Q G" 


^ T r ” T ■ 

-(k+l)  = -(k+l)  - (k+l)  L-(k+l)  - (k+l)  - (k+l) 


-(k+l)  - (k+l)  -(k  1) 


(9) 


^(k+1)  " - (k+1)  ■ -(k+1)  -(k+1)  ^ (k+1) 


Superscripts  (-)  and  (+)  denote  before  and  after  a measurement 
respectively.  ^ in  Equations  (5)  and  (6)  is  the  state  transition 
matrix  relating  the  states  at  time  instant  (k+1)  to  the  states 
at  time  instant  (k) . (See  page  22) . The  matrix  ^ in  Equation 
(6)  represents  the  strengths  of  the  white  noises  that  are  added 
to  the  last  (or  highest)  derivative  states.  (See  page  20). 

R in  Equation  (7)  is  a scalar  denoting  the  uncertainty  in  the 
DME  measurements.  (See  page  31)-  residual,  Ae  in 

Equation  (8) , is  the  difference  between  an  actual  DME  measurement 
and  what  the  estimator  predicts  it  to  be.  (See  page  27)  • The 


~(k+l)  Equations  (7)  and  (9)  is  the  partial  of  Equation 


(4)  with  respect  to  the  estimator  states,  evaluated  at  . 


(See  page  27)  • 

First,  Equations  (5)  and  (6)  are  used  to  propagate  the 
states  and  state  covariances  up  to,  but  not  including,  the  next 
measurement  as  and  respectively.  Second,  the 

Kalman  gain  matrix,  Z(k+1) ' calculated  from  the  known  state 
covariance  propagation.  It  should  be  noted  that,  unlike  the 
linear  Kalman  filter,  the  EKF  covariance  matrix  cannot  be 
precomputed  because  is  a matrix  of  partial  derivatives 

evaluated  at  , and  thus  requires  knowledge  of  the  measure- 

ment history.  Finally,  after  the  incorporation  of  each  measure- 

A 

ment,  the  estimator  updates  the  states  to  and  the 

state  covariances  to  via  the  Kalman  gain. 

System  Dynamics 

System.  The  system  in  this  study  is  an  aircraft  flying  at 


li< 


J 


some  known  altitude  and  known  initial  position  and  velocity. 

The  propagation  equations,  Equations  (5)  and  (6),  require  a 
knowledge  of  the  dynamics  of  this  system.  In  other  words, 
all  system  states  - position,  velocity,  and  so  forth  - are 
propagated  through  space  obeying  a set  of  dynamic  equations 
that  relates  the  system  to  its  environment. 

Approximations . To  avoid  a set  of  dynamic  equations  that 
is  too  complex  for  real-time  applications,  the  equations  are 
simplified  or  approximated  whenever  feasible.  One  simplification 
stems  from  the  fact  that  the  aircraft  flies  at  a known 
altitude.  The  availability  of  an  autopilot  with  an  accurate 
altimeter  enables  the  estimator  to  omit  the  vertical  states. 

Latham  has  shown  that  typical  altimeter  errors  degrade  the 
position  CEP  by  less  than  57o.  (Ref  5:152).  The  vertical  direction 
is  identified  with  the  up  direction  in  the  local  geodetic 
frame.  Since  the  vertical  state  is  absent,  the  estimator  assumes 
no  uncertainty  in  that  direction,  and  only  latitude  and  longitude 
states  are  estimated. 

In  the  9 - state  estimator,  an  assumption  that  the  system 
maintains  roughly  constant  rate  of  change  in  acceleration 
is  another  approximation.  Just  how  "rough"  depends  on  both 
pilot  control  and  natural  disturbances.  The  time  between 
measurements  (.05  to  .5  seconds)  is  short  enough  to  justify 
that  such  disturbances  cannot  significantly  alter  the  rate  of 
change  in  acceleration.  For  the  7 - state  estimator,  the  above 
approximation  is  moved  to  the  next  lower  derivative  level. 

Since  the  7 - state  estimator  does  not  model  jerk,  the 
acceleration  is  assumed  essentially  constant  between  measurements. 
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Another  simplification  for  the  9 - state  estimator  is 


that  each  jerk  state,  denoted  by  n(t),  is  modelled  as  the 
output  of  a first-order  lag  driven  by  a white  Gaussian  noise 
(totally  random  noise)  , a>(t)  . Diagram  (a)  depicts  the  relation- 
ship between  a)(t)  and  ri(t)  in  general  Laplace  form; 


(^(s) 


1 

s + a 


n(s) 


This  effectively  claims  that  disturbances  in  jerk  are 
exponentially  time-correlated.  Again,  each  acceleration  state 
in  the  7 - state  estimator  is  treated  in  the  same  manner 
as  the  jerk  states  in  the  9 - state  estimator. 

All  three  of  these  simplifications  are  employed  in  the 
estimator  design.  The  justification  of  such  simplifications 
is  found  via  estimator  performance. 

Estimator  States.  Before  the  set  of  equations  that  relate 
the  states  at  time  instant  (k+1)  to  the  states  at  time 
instant  (k)  can  be  designed,  the  various  states  of  the  estimator 
are  propagated  from  time  instant  (k)  to  (k+1)  with  coordinates 
of  latitude  and  longitude  in  units  of  radians.  Trajectory 
states  for  the  7 - state  estimator  are  chosen  as  position, 
velocity,  and  acceleration  in  both  latitude  and  longitude 
(i.e.  north  and  east)  directions.  The  9 - state  estimator  has 
the  jerk  state  added  in  both  directions.  Another  state  (of 
both  estimators)  is  the  DME  bias  which  is  also  propagated  and 
updated  in  each  measurement  interval.  If  latitude  and 
longitude  are  denoted  by  (J)  (subscripted  "g"  means  geodetic) 

O 

and  X,  respectively,  and  position,  velocity,  acceleration,  and 
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I jerk  are  denoted  by  P,  V,  A,  and  n , respectively,  the  states 

‘ of  the  two  estimators  can  be  summarized  as  in  Figure  2. 


Figure  2.  Estimator  States 


Dynamic  Equations.  If  the  states  in  Figure  2 are  represented 
in  an  array,  X,  the  differential  equations  relating  the  states 
to  each  other  can  be  written  as 


f- 


dX(t) 

dt 


= £ (X(t) , w(t)) 


(10) 


A discrete  equation  relating  the  values  of  to 

can  be  written  as 

^(k+1)  " - ^-(k)’  -(k)^ 

Estimator  states  at  time  instant  (k+1)  , or  , are 

a function  of  the  states  at  time  instant  (k)  , or  , 

and  a white  noise  vector  • — (k)  ^®P^®sents  the  deviations 

from  a constant  rate  of  change  in  acceleration  in  the  9 - state 
estimator  and  represents  the  deviations  from  a constant 
acceleration  in  the  case  of  the  7 - state  estimator. 

The  realtionship  between  a)(s)  and  ri(s)  in  each  direction 
for  the  9 - state  estimator  is  given  by  diagram  (a) . The 
value  of  "a"  in  this  diagram  is  set  to  zero  for  a pure  integrator. 
(It  should  be  noted  that  only  the  case  of  "a"  equal  to  zero  is 
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investigated  in  this  report,  but  other  values  of  "a"  are  suggested 
for  future  study).  With  "a"  equal  to  zero,  the  continuous 
stochastic  process  that  models  the  change  in  jerk  rate  as 
white  noise  is 

(12) 

where  n is  jerk  and  6j(t)  is  zero-mean  white  noise.  The 
discrete  model  for  the  expected  values  of  P,V,A,  and  n 
is  obtained  by  taking  the  expected  value  of  Equation  (12) : 


dn(t)  ^ Q 
dt 


(13) 


Given  the  best  estimates  of  the  states  at  t^j^^  ^^(k)  ^(k) 


^(k) , and  Equation  (13)  can  be  integrated  from  to 


n(t)  = n 


(k) 


(14) 


Integrating  Equation  (14)  from  t^^^  to  t yields 


A(t)  = 


(15) 


Expected  velocity  is  then  obtained  by  integrating  Equations 
(15): 


V(t)  +A^j^^(t  - ■ ^(k)^ 


(16) 


An  integration  of  Equation  (16)  provides  an  expression  for 
expected  position: 


" ^(k)  '^(k)^'^  " ^(k)'*’^^(k)^^'^(k)  ?’^(k)^^‘^(k)^^ 


(17) 
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Evaluating  Equations  (14) , (15) , (16) , and  (17)  at  time 
instant  (k+1)  uields 


/■ 


'’(k+1) 

'’(k) 

(18) 

^(k+1) 

A ^ 

" ^(k)  ’^(k)^^^^ 

(19) 

(k+1) 

A, 

" ’^(k) 

A 

* '^(k) 

A 

At  + % At^ 

(20) 

(k+1) 

A 

= ^k) 

A 

A A 

At  + % A(^^At^  + 1 At 

(21) 

where 

“ ^(k+l)‘^(k) 

Equations  (18)  through  (21)  are  the  mathematically  exact 

, and 


(22) 


integrations  of  Equations  (13)  from  to 


provide  expected  values  of  P,V,A,  and  n at  t 


(k+1) 


given  their 


best  estimates  at  Essentially,  the  mean  state,  values 

are  propagated  from  to  . 

Equations  (18)  through  (21)  can  be  expressed  as  the 
propagation  equations  between  measurements  for  the  9 - state 
estimator  when  the  expected  state  values  at  are 

recognized  as  the  estimator's  best  estimate  of  the  states  at 
t(k  • In  summary,  the  dynamic,  or  propagation,  equations 
of  the  9 - state  estimator  are 


■'+  ■'+  2 
P.  /I  , n \ = Pj.  /I  \ V.  V At  + % A.  /I  \ At  + 

<))g(k+l) 


O 


^<|.g(k+l)  “ %g(k)  %g(k) 


%g(k+l)  “ %g(k) 


^X(k+1)  " ^X(k)  ■*■  ^X(k)^^  ■’■  ^ ^X(k)^^^  ■•■  F ^X(k)^^^ 


^x(k+i)  " '^x(k)  ■*■  ^X(k)^^  ^ ’^X(k)^^^ 


^X(k+1)  ^X(k)  ’^X(k) 


^g(k+l)  ^({.g(k) 


^(k+1)  ’^X(k) 


To  model  the  uncertainty  in  the  propagation  equations, 
a first-order  approximation  of  Equation  (12)  yields 


^(k+1)  '^(k)  ‘"(k) 

which  can  be  expressed  in  latitude  and  longitude  directions: 


+ At  u) 


♦g(k) 


"X(k+l)  ' "xCk)  ■•■  "x(k) 

A similar  approach  is  used  for  the  7 - state  estimator 
except  the  noises  are  added  to  the  acceleration  states. 

The  propagation  equiations  for  the  7 - state  estimator  are  as 
follows : 


1 ' 


ry 
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^<{.g(k+l)  = ^^g(k)  ^ ^^g(k) 


^<^g(k+l)  ~ '^lg(k)  ^^g(k) 


%g(k+l)  ~ %g(k) 


X(k+1)  '^A(k)  ^ ''x(k) 


P^/1.l^\  Pi/i.N  V-v/i.\  At  + % A^  At 


^X(k+1)  '^X(k)  ■'■  ^X(k) 


\(k+l)  - ^X(k) 


(34) 


(35) 


(36) 


(37) 


(38) 


(39) 


Entering  the  noise  at  the  acceleration  level  and, 
again,  using  a first-order  approximation  provides  expressions 
for  the  uncertainty  in  the  propagation  equations: 


%g(k+l)  %g(k)  + ^%g(k) 


^X(k+1)  ■ ‘^X(k)  ‘^X(k) 


(40) 


(41) 


An  additional  state  is  augmented  to  the  estimator  states 
(of  both  estimators)  due  to  the  presence  of  the  DME  range 
measurement  bias.  Measurement  bias  is  modelled  as  an  integrator 
with  no  white  noise  input  and  a random  initial  condition. 
Integrating  a random  initial  condition  produces  what  is  known 
as  a random  constant  or  bias,  b^(t). 
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a- 


h' 


Initial 

Condition 


dt 


b.(t) 


(b) 


The  bias  differential  equation  from  diagram  (b)  is 


db,(t) 


= 0 


dt 


(42) 


Integrating  Equation  (42)  from  to  yields 


^i(k+l)  ^i(k) 


(43) 


The  subscript  (i)  denotes  a bias  associated  with  station  (i) . 

In  changing  from  station  to  station,  the  initial  condition  in 
diagram  (b)  must  be  reset  to  correspond  to  the  correct  station. 
The  method  of  separating  the  bias  according  to  station  is 
covered  in  the  section  on  "Bias  Calculation"  in  this  chapter. 

A Gaussian  random  variable  model  is  used  for  the  initial 
condition  with  zero  mean  and  a variance  consistent  with 


typical  bias  errors  given  in  Chapter  II.  In  this  study, 
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bias  variance  is  set  at  10  rad  , or  432  ft  . The  use  of 


different  bias  variances  are  not  attempted  in  this  report 
but  are  suggested  for  further  study. 

Covariance  Propagation 

The  state  propagation  of  the  previous  section  provides 
the  estimator  with  a state  transition  matrix,  ^ is  the 
matrix  representation  of  the  dynamic  equations,  and,  for  the 
9 - state  estimation,  $ is  denoted  by  the  following  9x9 
matrix: 


2? 


J 


At  /, 


3 

Af^/, 


At  /. 


At  At'^/2  0 Af^/^ 

1 At  0 At^/, 


^ for  the  7 - state  estimator  is  obtained  in  the  same  manner. 
The  white  noise  coefficients  from  Equations  (32)  and  (33)  are 
represented  in  a separate  matrix,  G; 


00000  At  00 
0 0 0 0 0 0 At  0 


The  covariance  is  propagated  between  measurements  as 

^ A 

r(k+i)  - i p^k)  i''  + s 9 e ^ (6) 

/s 

The  diagonal  terms  of  the  9x9  covariance  matrix,  » 

give  variances  for  state  errors  before  each  measurement.  The 
off-diagonal  terms  yield  an  estimate  of  the  correlation  between 
each  state.  Near  perfect  initial  conditions  (position  and 
velocity)  are  assumed  for  this  estimator  so  is  initiated 
with  small  starting  values. 

The  matrix  ^ in  Equation  (6)  is  a 2 x 2 matrix  denoting  the 
strengths  of  the  white  noises  o)  and 
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S = 


(46) 


The  values  of  the  elements  in  equation  (46)  are  fixed  to 
reflect  disturbances  in  jerk  for  the  7 - state  estimator  and 
disturbances  in  jerk  rate  for  the  9 - state  estimator. 

The  state  propagation,  Equation  (5),  and  state 
covariance  propagation,  Equation  (6) , are  used  to  propagate 
the  states  and  covariances  between  each  measurement.  Now 
the  estimator  is  ready  to  use  incoming  measurements 
to  update  the  states  and  covariances  at  each  measurement 
instant.  Before  the  update  equations  are  employed,  the 
estimator  must  be  able  to  predict  the  value  of  each 
measurement.  Then  a differencing  of  the  prediction  and  the 
actual  measurement  (i.e,  the  residual)  can  be  used  as  a 
means  to  update  the  states  and  covariances. 

Measurement  Prediction 

Introduction.  Once  the  state  values  have  been  propagated  up  to 
the  end  of  the  measurement  interval,  the  estimator  predicts  a 
value  of  the  next  measurement.  The  prediction  is  the  "best 
guess"  at  what  the  measurement  should  be  and  is  based  on  all 
prior  knowledge  of  the  states  up  to,  but  not  including  the  next 
measurement.  DME  station  coordinates  and  the  current  estimate 
of  aircraft  coordinates  are  assigned  to  a common  Cartesian  coordinate 
frame.  The  distance  (or  measurement)  equation.  Equation  (4),  is 
then  employed  to  yield  a value  for  the  measurement  prediction. 


Essentially,  the  measurement  prediction  is  divided  into 
three  parts:  DME  station  information,  transformation  of 
latitude,  longitude,  and  height  to  a Cartesian  coordinate 
frame,  and  finally  the  prediction  of  the  measurement  value 
itself. 

DME  Station  Information.  Incorporation  of  the  channel 
information  is  accomplished  with  a subroutine,  containing 
all  the  latitudes,  longitudes,  and  heights  of  all  the  DME 
stations  the  aircraft  might  use  on  a particular  flight.  A 
simple  "table-lookup"  routine  is  used  for  associating  the 
correct  channel  number  to  the  respective  station  coordinates. 

This  station  information  is  obtained  from  DoD  Flight  Information 
Publication  (IFR-Supplement , issued  every  eight  weeks). 

Latitude  and  longitude  are  in  units  of  degrees  in  the  supplement 
and  must  be  converted  to  radians  for  the  estimator.  The 
latitudes  and  longitudes  in  the  supplement  correspond 
to  those  on  all  Air  Force  maps  (geodetic  latitude) . DME 
station  information  is  inserted  into  the  subroutine  prior  to 
each  flight.  The  estimator  then  extracts  this  information 
from  the  subroutine  in  real-time. 

Coordinate  Transformation.  In  order  to  calculate 
a predicted  range  value,  the  latitudes,  longitudes,  and 
heights  of  the  DME  station  and  aircraft  must  be  transformed 
to  a common  Cartesian  coordinate  frame.  Range  can  then  be 
computed  using  the  geometric  distance  equation.  Equation  (4),  for 
computing  distances  between  two  points  in  three  dimensional 
space.  As  mentioned  before,  an  appropriate  Cartesian  frame 
for  the  estimator  is  an  XYZ  frame  centered  at  the  earth's 


'i 

i 


i 

i 


25 


center  with  the  Z-axis  through  the  north  pole,  the  Y-axis 
through  the  Greenwich  Meridian  and  equator,  and  the  X-axis 
forming  a right-handed  coordinate  system.  Formulas  that 
express  X,  Y,  and  Z Cartesian  coordinates  as  functions  of 
geodetic  latitude,  longitude,  and  height  can  be  written  as 

X = ^1 

g 

Y = f2  (P^  . P;^,  Hgt)  (47: 

Z = f3  (P^  , P^,  Hgt) 

where  Hgt  is  height.  Appendix  A presents  this  set  of  equations 
in  detail.  These  equations  are  used  to  obtain  X,Y,  and  Z 
coordinates  for  both  station  and  aircraft  locations. 

Measurement  Prediction.  A knowledge  of  the  station  and 
aircraft  coordinates  now  allows  the  range  estimate  to  be 
calculated  from  the  measurement  equation.  Equation  (4) , using 
the  zero  mean  value  of  the  white  noise,  v. 


e = (X  -X  )^+  (Y  -Y  )^  + (Z  -Z  ) h + b. 
L'-  u s'  ^ u s ^ u s'  J 1 


(k+1) 


e represents  the  predicted  range  value;  X , Y , and  Z 

s s s 

^ Ak  y\ 

are  the  station  coordinates,  and  X^,  Y^,  and  Z^  are  the 
estimated  aircraft  coordinates.  An  estimated  bias, 

is  added  directly  to  the  measurement  prediction  equation. 

A 

A comparison  of  the  predicted  measurement,  e,  with  the 
actual  incoming  DME  measurement  provides  the  estimator  with 
a residual  for  updating  the  states  and  covariances  for 
each  measurement.  Before  proceeding  into  the  updating  process. 
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The  cutoff  point  has  both  a lower  and  an  upper  limit.  The 
residual  cutoff  point  is  made  small  enough  to  ignore  bad 
measurements,  but  it  must  also  be  set  large  enough  to  allow 
for  aircraft  trajectory  changes.  (Trajectory  changes  such 
as  180“  turns  will  tend  to  temporarily  increase  residual 
size  (See  Chapter  V)). 

The  updating  process  is,  in  essence,  a method  for 
changing  the  propagated  states  by  amounts  proportional  to 
the  residual.  Obtaining  an  exact  proportion,  or  Kalman  gain, 
constitutes  a central  issue  of  Kalman  filtering,  and  requires 
a linearized  measurement  equation  to  generate  the  desired  H 
matrix  for  Equations  (7)  and  (9) . Since  Equation  (4)  is 
nonlinear,  some  method  for  approximating  this  equation  as 
linear  constitutes  the  next  step  in  the  estimation  process. 

Linearized  Measurement  Equation 

The  residual  computation  of  the  last  section  employs  the 
exact  non-linear  measurement  equation,  but,  to  permit  matrix 
operations  in  the  update  equations,  the  measurement  equation 
must  be  linearized.  This  linearization  allows  the  measurement 
equation  to  be  represented  by  a linear  measurement  matrix, 
or  H matrix.  The  H matrix  forms  the  crux  for  the  update 
Equations  (7),  (8)  and  (9). 

To  calcualte  the  H matrix,  the  complete  expression  for  e, 

" - [«u  - + V (50) 

can  be  expressed  as  a first  - order  Taylor  series  linearized 
about  the  best  estimates  of  the  states.  As  long  as  the  estimator 
remains  close  to  the  true  state  values,  H will  adequately 
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represent  the  non-linear  measurement  equation,  A first- 
order  Taylor  series  of  Equation  (50)  yields  e as  a linear 


function  of  X^,  Y^,  Z^,  X^,  Y^,  Z^,  and  v: 


y\  ^ A /V 


e(X,Y,Z,X,Y,Z,b.,v)=e(X,Y.Z,X,Y, 
' u’  u’  u’  s’  s’  s’  i’  ^ u’  u’  u s’  s’ 


2s.  b.) 


(^u  -^u>  W 


(Y  - Y ) 
' u u 


Yu. 


, ^ 

32. 


Z ) + 

V 3b. 


(b . - b . ) + V 
' 1 1.' 


X .Y  ,Z 
u u u 


Errors  in  station  coordinates,  X , Y , and  Z , are  assumed 

s s s’ 

to  be  negligible. 

Equation  (51)  can  be  written  equivalently  as  a 
function  of  P.  , P, , height,  b.,  and  v using  Equation  (47); 

<Py  A 1 


e(P^  . P;^,  Hgt,  b.,  X^.Y^,  Z^)  = e (P^  .P^.Hgt,  b^.X^.Y^.  Z^) 

s s 


({•g  g g ^ 


-(k+1) 


-(k+1) 


If.,  (^i  - bi)  + V 


-(k+1) 


r 


y; 


K 


where  errors  in  height  are  assumed  negligible  and 

is  the  best  estimate  of  the  states. 

An  investigation  of  Eq  (52)  shows  that  the  equation  is 
in  the  form 

e - e = Ae  = [H]  [X  - X]  + v 
where  H is  given  by 


3e 


H = 


3P, 


3e 

3e 

o 

o 

o 

o 

o 

o 

3PA 

3b^ 

(53) 


-(k+1) 


-(k+l) 


Kk+1) 


The  units  of  the  H terms  are  feet/radian  since  the  residual, 
Ae,  is  in  feet  and  the  states  are  in  radians. 

The  actual  evaluation  of  the  partials  of  Equation  (53) 
involves  the  substitution  of  Equation  (47)  into  the  measure- 
ment prediction  Equation  (50) . The  details  of  this  operation 
are  presented  in  Appendix  B. 

The  bias  is  added  directly  in  Equation  (50) , and 
consequently  an  incremental  change  in  bias  produces  the 
same  incremental  change  in  range.  In  this  light,  it  is 

evident  that  3e/3b.  = 1 when  e and  the  bias  are  in 

1 

common  units.  However,  since  e is  in  units  of  feet  and 
bias  in  units  of  radians,  3e/3b^  is  in  units  of  feet/radian. 
Because  longitude  lines  converge  as  the  latitude  angle 
increases,  the  number  of  feet  per  radian  cannot  be  calculated 
using  the  conversion  factor  of  60  nautical  miles  per  degree. 
Both  latitude  and  longitude  components  are  involved  in  the 
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f: 


direction  of  the  range,  thus  making  each  degree  contain  less 
than  60  nautical  miles  in  the  calculation  of  3e/9b^. 

The  orthogonality  of  the  directions  from  which 
latitude  and  longitude  are  measured  yields  a simple 
method  for  obtaining  the  value  of  the  bias  partial.  Figure  3 
shows  that  3e/3b^  contains  components  of  3e/3P^  and  3e/3P^. 
Applying  the  Pythagorean  theorem  provides  an  expression  for 


3e/3b^; 


3e 

" L 37. 

X ({) 


3e 


2 ^ 


] 


(54) 


g 


The  dependency  of  3e/3b^:  on  3e/P^ 


and  3e/P,  leaves 


8 


the  estimator  only  the  latter  two  partials  to  determine 


LATITUDE 


97, 


Figure  3.  Dependence  of  3e/3bj^  on  and  3e/3P^  in 


H Matrix. 
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The  equations  needed  to  evaluate  9e/8p,  and  9e/3P, 
always  depend  on  the  current  best  state  estimates- a property 
of  any  extended  Kalman  filter.  Essentially,  for  each  new 
measurement  a new  set  of  nominal  conditions  is  used  to 
solve  the  first-order  Taylor  series.  The  H matrix  will 
therefore  change  its  values  from  one  measurement  to  the  next. 
For  each  measurement,  the  algorithm  recalculates  an  H and 
utilizes  it  in  the  state  and  state  covariance  update  process. 
Kalman  Gain  and  Update  Equations 

The  update  process  consists  of  the  implemention 
of  Equations  (7),  (8),  and  (9).  The  update  portion  is 
accomplishable  only  after  the  states  and  state  covariances 
are  propagated  and  after  the  residuals  and  linearized 
measurement  equation  are  obtained.  The  algorithm  divides 
the  update  process  into  three  separate  steps:  the  Kalman 
gain,  the  state  update,  and  the  covariance  update. 

The  proportion  by  which  the  propagated  estimates  are 
changed  after  the  incorporation  of  each  measurement  is 
called  the  Kalman  gain.  Equation  (7)  is  repeated  here 
for  convenience : 


^(k+1)  = -(k+1)  -(k+1) 


[ 


^(k+D-  (k+1) -(k+1) 


(7) 


The  Kalman  gain,  , essentially  determines  how  much 

"faith"  the  estimator  has  in  the  measurements  with  respect 
to  the  state  propagations.  Since  changes  from 

measurement  to  measurement,  also  changes. 

The  value  of  R represents  the  uncertainty  of  incoming 
DME  measurements;  a low  R value  reflects  accurate  measurements. 
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The  DME  error  model  presented  in  Chapter  II  gives  a "ball- 
park" range  for  R.  Proper  R values  are  best  obtained  by  a 

process  called  filter  tuning.  However,  in  this  report, 

2 

R is  set  at  1000  feet  . So  allow  for  changes  in  R,  R can 
be  made  a function  of  time,  represented  by  R(t) . It  may  be 
advantageous  to  have  a changing  R to  allow  for  time-varying 
noise  characteristics  in  the  real  world  system  (such  as 
due  to  changing  ranges  from  DME  stations) , 

Sometimes  the  Kalman  gain  is  referred  to  as  a 
weighting  matrix  because,  in  a sense,  the  algorithm 
"weighs"  the  residuals  with  the  Kalman  gain.  The  weighted 
residuals  are  then  added  to  the  best  prior  state  estimates 
to  achieve  a new  updated  estimate  of  the  states.  This  is 
accomplished  by  Equation  (8) ; 


-(k+1)  -(k+1)  + -(k+1) 


(8) 


where  represents  the  updated  estimates  and  Ae  is  the 

residual . 

The  propagated  state  covariance  is  updated  in  the  same 
fashion  using  Equation  (9) : 

^ I tm  /N 

-(k+1)  = - (k+1)  ■ -(k+1)  -(k+1)  - (k+1)  (9) 


As  in  the  state  case,  a new  updated  estimate  of  the  state 
covariance  matrix  is  formed. 

Bias  Calculation 


The  bias  calculation  is  an  example  of  an  ad  hoc  procedure 
to  be  incorporated  in  the  estimator.  The  bias  is  composed 
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of  random  constant  errors  in  both  the  airborne  equipment  and 
the  DME  stations.  Because  different  stations  are  acquired, 
the  bias  will  change  from  measurement  to  measurement. 

However,  the  bias  must  remain  unique  to  each  corresponding 
station.  Some  type  of  separation  and  bias  variance  reinitialization 
method  is  needed  to  maintain  properly  corresponding  station 
biases . 

This  need  for  the  algorithm  to  separate  bias  according 
to  station  has  led  to  the  following  method.  The  bias  is 
stored  in  an  n x 1 array,  where  n is  equal  to  or  greater 
than  the  number  of  stations  acquired  during  the  flight. 

Bias  values  obtained  from  the  state  update  equations 
are  stored  in  different  elements  of  the  array  for  different 

channels.  If  a channel  is  acquired  more  than  once  (which 

is  necessary  for  a good  bias  estimation) , then  only  the  bias 
array  element  corresponding  to  that  station  is  called  on 
as  the  best  prior  estimate  of  the  bias  state.  Using 
mmm  as  an  index  and  BI  as  the  bias  array,  the  estimator 
stores  the  bias  as 

BI  (mmm) 

BI  (1)  ->■  mmm  = 1 for  Channel  A 

BI  (2)  ->■  mmm  = 2 for  Channel  B 

BI  (3)  ^ mmm  = 3 for  Channel  C 

etc . 

A,  B,  and  C denote  arbitrary  channels.  BI  is  dimensioned 
at  least  the  maximum  amount  of  stations  acquired  in  the 
flight.  The  estimator  returns  to  the  appropriate  mmm  value 
if  that  corresponding  station  is  called  on  again,  or  if, 
in  the  case  of  a new  station,  mmm  is  incremented  by  one. 
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For  example,  after  the  estimator  accepts  40  measurements  and 
Channel  B is  called  on  again  for  a range  measurement,  the 
BI  array  element  associated  with  mmm  = 2 is  updated  and 
stored  again  in  BI  (2) . 

The  separation  technique  must  also  include  a method  for 
reinitializing  bias  variance  for  each  time  a new  DME  station 
is  acquired.  Each  time  a new  station  (not  previously  used) 
is  acquired,  bias  variance  is  reinitialized  to  the  maximum 
value.  If  a station  has  been  previously  acuqired,  the  bias 
variance  is  reset  to  the  latest  value  corresponding  to 
that  station.  This  effectively  claims  that  a particular 
station's  bias  uncertainty  is  greatest  for  the  first  acquisition 
of  that  station.  As  that  station  is  used  again  and  again, 
bias  uncertainty  will  generally  decrease.  (See  Figure  15 
/;  in  Chapter  IV)  . The  actual  program  in  Appendix  C shows  the 
detailed  procedure  for  treating  the  bias  separation. 

Truth  Model  Data 

Data  from  an  aircraft  trajectory  is  necessary  to  test 
and  analyze  the  estimator  design.  Measurements  from  a 
simulated  flight  or  measurements  recorded  from  an  actual 
flight  can  be  used  to  test  the  estimator.  In  the  case  of 
a flight  simulation,  the  measurements  to  DME  stations  must 
be  corrupted  with  a noise  generator  (random  number  generator) . 
However,  recorded  measurements  taken  from  an  actual  flight 
provide  the  best  way  to  test  the  estimator. 

Raw  DME  measurement  data,  CIRIS  output  data,  and 
FASTMAP  filter  data  for  a test  flight  were  stored  on  tapes. 
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I Channel  information,  DME  measurements,  and  system  time  (IRIG  time) 

I from  the  FASTMAP  filter  tape  are  used  as  the  input  for  the 

I estimator.  The  use  of  FASTMAP  tape  allows  easier  comparison 

f of  FASTMAP  results  with  the  real-time  design  estimator. 

(See  Chapter  IV).  Essentially,  the  range  data  on  the 
FASTMAP  tape  serves  as  a replacement  for  simulated  system 
data. 

The  geometry  of  the  station  locations  relative  to  the 
aircraft  is  important  to  the  accuracy  of  the  estimated 
aircraft  position.  A more  spread  out  distribution  of  stations 
is  preferred  over  a situation  where  all  the  available 
stations  are  in  one  direction.  (See  Figure  4 and  Figure  5). 
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Figure  4.  Relative  Station  and  Aircraft  Locations 
Stations  in  Same  General  Direction. 
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In  Figure  4,  the  stations  are  nearly  colinear,  and  normal 
to  the  aircraft  trajectory.  Therefore  position  accuracy 
along  the  trajectory  is  much  worse  than  position  accuracy 
normal  to  the  trajectory.  This  results  in  a locus  of  constant 
likelihood  in  the  shape  of  a highly  eccentric  ellipse;  constant 
likelihood  lines  are  represented  by  the  dotted  lines. 


Figure  5.  Relative  Station  and  Aircraft  Locations: 
Stations  Radially  Distributed. 


When  the  distribution  of  stations  is  spread  out  as  in 
Figure  5,  errors  are  more  radially  distributed.  As  the 
flight  progresses  and  the  number  of  station  acquisitions 
increases,  the  probability  of  better  station  distributions 
increases.  The  first  30  minutes  of  the  FASTMAP  filter  data 
shows  the  use  of  20  different  DME  stations.  Figure  6 
illustrates  the  distribution  of  these  20  stations  relative 
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to  the  aircraft.  The  number  labels  in  the  figure  correspond 
to  the  channels  of  the  DME  stations.  "A"  is  the  start  of 
the  flight  and  "B"  is  the  aircraft  location  after  30  minutes 
of  flight. 


The  Overall  Estimator 


The  estimator  is  designed  to  produce  the  best  estimate 
of  the  states  by  first  propagating  the  states  from  time 
instant  (k)  to  (k+1) . The  state  update  is  accomplished  with 
a new  measurement  at  (k+1) . (k)  is  then  incremented 

by  one  and  the  procedure  is  repeated.  The  next  propagation 
begins  with  the  last  interval's  updated  estimates,  and  a 
new  measurement  updates  the  estimates  again.  In  transferring 

from  one  time  interval  to  the  next,  the  estimator  renames 

''+  ■'+ 
the  states,  , and  covariances,  . with  and 


''+ 

5 (k) 

(k» 


respectively  (in  other  words,  the  process  iterates  on 
This  process  is  represented  by 


^(k)  ^ ^(k+1) 

A 

^Ik)  i(k+i) 


The  overall  estimator  flow  chart  for  the  7 - state  estimator 
is  shown  in  Figure  7 . 
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PROPAGATE  COVARIANCE  MATRIX: 


Take 


Kalman  gain  and  update 


Figure  7.  7 - State  Estimator  Flow  Chart. 


IV . Estimator  Performance  Compared  to  the  FASTI'IAP  Filter 


Initial  Conditions 

The  estimator  algorithm  employing  the  design  of  the 
previous  chapter  is  tested  by  using  range  and  channel  data 
from  the  FASTMAP  tape.  Trajectory  data  also  recorded  on  the 
FASTMAP  tape  are  the  trajectory  estimates  from  the  post- 
flight FASTMAP  filter.  A direct  comparison  between  the 
estimator  and  the  FASTMAP  filter  trajectories  provides  a 
useful  evaluation  of  estimator  performance.  Initial  position 
and  velocity  for  the  estimator  are  obtained  from  the  initial 
position  and  velocity  of  the  FASTMAP  filter.  The  test  flight 
available  from  the  FASTMAP  test  begins  at  latitude  .606342425 
NORTH  (radians)  and  longitude  1.845317956  l-JEST  (radians).  A 
constant  height  of  31524  feet  is  maintained  by  the  aircraft 
autopilot.  The  aircraft  is  headed  in  a northerly  direction  at 
approximately  300  nautical  miles  per  hour.  Two  hundred  miles 
north  of  the  initialization  of  the  FASTMAP  test,  the  aircraft 
makes  a 180“  heading  change  and  returns  to  the  initial 
latitude  and  longitude.  Initial  conditions  are  as  follows: 


^<l>^(o)  ” .606342425  radian  = 34.7409  degrees  (55) 
V^^(o)  = 2.2666  X 10  ^ rad/sec  = 279.06  nm/hr  North  (56) 
A^^(o)  = 0.0  rad/sec^  = 0 ft/sec^  (57) 


P 


+ 

Uo) 


1.845317956  radian  = 105.7289  degrees 


(58) 


i 


+ 

^(O) 

= 9.2103  X 10  ^ rad/sec  = 90.72  nm/hr 

West  (59) 

A + 
^(o) 

= 0.0  rad/sec^  = 0 ft/sec^ 

(60) 

^i(o) 

= 0.0  radian  = 0 feet 

(61) 

9-state 

estimator  initializes  the  jerk  states 

at  zero: 

= 0.0  rad/sec^  = 0 ft/sec^ 

(62) 

A 

^X(o) 

=0.0  rad/sec^  = 0 ft/sec^ 

(63) 

Initial  covariance  values  depend  on  the  uncertainty  of 

the  initial  values  of  the  FASTMAP  filter.  Low  uncertainty 

for  these  initial  conditions  is  assumed;  therefore,  the 

diagonal  terms  of  the  initial  covariance  matrix,  or  ^ , 

are  set  at  small  values.  In  other  words,  initial  covariance 

is  made  somewhat  smaller  than  the  anticipated  steady  state 

covariances  for  the  estimator.  Position  variances  are 

chosen  as  10  radians  squared,  or  4.3  feet  squared,  which  is 

less  than  the  anticipated  steady  state  estimator  covariances 

(on  the  order  of  10^  feet  squared,  consistent  with  the 

estimator  goals).  Velocity  and  acceleration  initial  variances 

are  assigned  in  the  same  manner.  The  correlations  (off- 

diagonal  terms)  between  the  states  are  initially  assumed  to  be 

zero.  After  the  estimator  begins,  both  the  diagonal  terms  and 

off-diagonal  terms  will  increase  except  for  the  bias  variance. 

The  bias  variance,  ^ (9,9)  in  Equation  (64),  is  initially 
-12 

set  at  10  radians  squared  or  432  feet  squared,  a value  more 
or  less  consistent  with  DME  range  bias  errors  as  outlined  in 


in  Chapter  II,  After  the  estimator  begins,  bias  variance  in 

the  estimator  is  expected  to  decrease  due  to  incorporating 

measurement  information. 

The  P matrix  for  the  9-state  estimator  is 
— o 


lo  = 


10"^^(rad)^ 


10~^^(rad/sec) ^ 

10^^(rad/sec^)^ 
10^^ (rad) ^ 


10  ^^(rad/sec)^ 

10  (rad/sec^) ^ 

10  (rad/sec^) ^ 


10"^^  (rad/sec-") 


3^2 


10‘^^(rad)2 


(64) 


The  P^  matrix  for  the  7-state  estimator  is 


4 = 


lO'^^(rad)^ 


10’^^(rad/sec)^ 


10 


■15 


2^2 


(rad/sec  ) 
10"^^(rad)^ 


10  (rad/sec) ^ 

10"^^(rad/sec^)^ 

10"’-^(rad)^ 


(65) 
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The  values  of  the  elements  of  Equation  (6)  are 
fixed  to  reflect  disturbances  in  jerk  for  the  7-state 
estimator  and  disturbances  in  jerk  rate  for  the  9-state 
estimator.  During  the  testing  of  the  estimator  (Chapters 
IV  and  V),  such  disturbances  are  assumed  to  be  small;  thus 
low  values  for  the  ^ elements  are  used.  The  units  of  are 
obtained  from  the  equation 


E {w(t)  u'^Ct  + t)}  = ^(t)  5(t)  (66) 

For  the  7 -state  estimator,  w is  added  directly  to  acceleration 

3 

through  an  intergrator  and  has  units  of  rad/sec  . The  units 
of  (the  5 matrix  for  the  7-state  estimator)  are  obtained 
from  a units  breakdown  of  Equation  (66) : 

(rad/sec^) (rad/sec^)  = ^(t)l/sec  (67) 

or  is  in  units  of  radians  2/seconds^.  A similar  approach 

is  employed  for  the  9-state  estimator  to  yield  the  units  of 
2 7 

as  rad  /sec  . For  this  study  (but  not  shown  here),  several 

different  values  of  Q were  initially  attempted,  and  the  best 

2 5 

5 values  seemed  to  be  those  smaller  than  1 ft  /sec  for 
2 7 

and  1 ft  /sec  for  Nevertheless,  it  should  be  emphasized 

that  fine  tuning  is  not  attempted  in  this  study,  and,  there- 
fore, only  "ball  park"  figures  for  3.  suggested.  A 

2 2 

conversion  from  ft  to  rad  yields 

diag  < l(ft^/sec^) (1/6076  nm/ft)^(l/60  deg/nm)^ 

(1/57  rad/deg)^  = 2.3  x 10"^^  rad^/sec^ 


i'lii'ir'h  I 


1 

I 


i 
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diag  {^g}  < 2,3  X 10  rad^/sec^ 


(68) 


The  ^ diagonal  values  for  both  estimators  are  chosen  as  lO”^^ 

2 5 -17  2 7 

rad  /sec  or  10  rad  /sec  and  are  stored  in  the  2x2 
matrix  as  follows: 

= S9  = 

1 
I 

Note  that  these  Q values  are  design  parameters  and  can  be 
changed , 


10 

0 


-17 


10 


-17 


(69) 


Comparison  of  7-  and  9-state  Estimators 

Both  the  7-state  and  the  9-state  estimator  positions 
and  velocities  are  compared  to  the  FASTMAP  filter  positions 
and  velocities  using  the  specified  initial  conditions.  Since 
the  actual  test  flight  initially  begins  straight,  level,  and 
with  zero  acceleration,  proper  estimation  results  should 
indicate  this.  The  actual  results  of  testing  the  estimators 
show  that  the  9-state  estimator  velocity  wildly  oscillates 
about  the  FASTMAP  velocity.  On  the  other  hand,  the  7 -state 
estimator  has  the  opposite  effect;  in  other  words,  its 
velocity  curve  is  even  smoother  than  the  FASTMAP  velocity  curve. 
(Note  that  R in  either  case  is  kept  constant  at  1000  ft  so 
that  effects  of  measurement  uncertainty  are  the  same.)  This 
phenomenon  eventually  leads  to  the  exclusion  of  the  9-state 
estimator  from  the  design  in  favor  of  the  7-state  estimator. 
However,  better  performance  may  be  attainable  for  the  9- state 
estimator  through  tuning.  In  fact,  in  this  light  no 
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conclusive  thoughts  can  be  drawn  from  this  comparison  except 
that  utilizing  the  current  values  of  ^ and  R,  the  7 -state 
estimator  outperforms  the  9-state  estimator. 

Estimator  positions  are  essentially  the  same  for  both 
the  7-  and  9-state  estimators.  For  example,  at  the  conclusion 
of  the  first  record  of  data  (42  DME  measurements  correspond 
to  one  record  of  data  and  the  entire  flight  is  represented  by 
over  500  records  on  the  FASTMAP  flight  data  tape) , updated 
position  state  values  of  both  estimators  compare  well  as  shown 
in  Table  II.  For  eay  comparison  in  Table  II,  FASTMAP  filter 
position  is  subtracted  from  the  positions  of  the  two  estimators. 
As  shown  in  the  table,  the  differences  in  both  cases  are  almost 
the  same,  thereby  proving  good  position  agreement  with  th. 

7-  and  9-state  estimators. 

Table  II.  Comparison  of  7-  and  9-State  Estimator  Positions 


-^Reference 

FASTMAP 

7-State 

9-State 

7-S  EST. 

9-S  EST 

Position 

Estimator 

Estimator 

Minus 

Minus 

(Radians) 

Position 

Position 

FASTl^P 

FASTMAP 

Direction^  . 

(Radians) 

(Radians) 

(Feet) 

(Feet) 

Latitude 

.606483606 

.606482181 

.606482155 

-29.6 

-30.0 

Longitude  1 

.845375889 

1.845366924  1.845366823 

149.0 

150.0 

When  velocities  of  the  9-state  estimator  and  FASTMAP 
are  compared,  an  interesting  result  is  obtained.  The  estimator 
velocity  states  oscillate  about  the  FASTMAP  filter's  velocity. 
The  magnitude  of  the  oscillations  are  quite  unacceptable 
ranging  from  -15  to  +15  nautical  miles  per  hour  from  the 


r 


iJ-' 


average  in  a 2 1/2  second  period.  To  illustrate  these 
oscillations,  longitude  velocities  of  the  estimator  and 
the  filter  are  compared  in  the  fifth  record  of  data.  (Using 
the  fifth  data  record  of  the  tape  should  insure  that  effects 
of  perfect  initial  agreement  are  gone  since  estimator  initial 
conditions  were  taken  from  the  FASTMAP  filter's  initial 
conditions.  Otherwise,  using  longitude  velocity  in  the  fifth 
record  is  an  arbitrary  choice.)  Table  III  tabulates  the 
time  - average  longitude  velocity  in  the  fifth  record  given 
by 

+ 

(70) 


42 

E 

i=l 


^(i) 
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and  the  high  and  low  longitude  velocities  in  the  same  interval 
for  both  the  9-state  estimator  and  the  FASTMAP  filter. 


Table  III.  9-State  Estimator  and  FASTMAP  Filter  Longitude 

Velocity  in  Fifth  Record 


>>,.^^anity 

Reference^ 

Time-Avg.  Long. 
Velocity 
(nm/hr) 

High  Long. 
Velocity 
(nm/hr) 

Low  Long . 

Velocity 

(nm/hr) 

9-State 

85.74 

125.83 

71.65 

Estimator 

FASTMAP 

84.74 

104.27 

71.68 

Filter 

The  table  shows  that  velocity  drops  for  this  particular 
period  of  time  to  equal  lows  for  both  the  9 -state  estimator 
and  the  FASTMAP  filter,  but  the  high  velocities  differ  by  over 
16  nautical  miles  per  hour.  Since  FASTMAP  tests  velocity 
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errors  are  +9.2  feet/second  CEP,  or  + 5.45  nautical  miles/ 
hour  CEP,  such  deviations  in  the  estimator  velocity  states 
as  illustrated  in  Table  III  cannot  be  accepted.  Therefore, 
the  performance  of  the  estimator  with  jerk  states  reveals 
that  the  9-state  estimator  must  be  rejected  (that  is,  with 
the  current  values  of  R and  . 

Table  IV  shows  a comparison  between  the  7-state 
estimator  and  the  FASTllAP  filter.  As  in  the  previous  test, 
the  comparison  is  accomplished  using  longitude  velocities  in 
the  fifth  record.  The  point  to  be  drawn  from  Table  IV  is  that 
the  estimator  velocity  varies  even  less  than  the  FASTMAP 
filter  velocity.  This  result  is  far  more  advantageous  than 
the  9-state  estimator  test  results  because  the  7-state 
estimator  velocity  better  reflects  the  fact  that  the  flight  is 
relatively  straight,  level,  and  has  near  zero  acceleration. 
Again,  it  must  be  emphasized  that  these  observations  are 
also  a function  of  tuning.  A suggestion  for  tuning  and 
testing  the  estimator  utilizing  different  trajectories  in  a 
Monte  Carlo  analysis  is  reserved  for  the  recommendations . 


Table  IV.  7- 

■State  Estimator  and 

Velocity  in  Fifth 

FASTMAP  Filter 

Record 

Longitude 

■“^^...Ouanity 

Time-Avg.  Long. 

High  Long . 

Low  Long . 

Velocity 

Velocity 

Velocity 

Reference 

(nm/hr) 

(nm/hr) 

(nm/hr) 

7-State 

Estimator 

84.52 

106.08 

84.25 

FASTMAP 

Filter 

84.74 

109.27 

71.68 
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Figure  8 shows  the  plot  of  longitude  velocity  for  the 
7-state  estimator  versus  the  9-state  estimator  for  the  fifth 
data  record.  The  large  magnitude  of  velocity  deviations  in 
the  9-state  estimator  relative  to  the  7-state  estimator  are 
apparent.  Due  to  these  velocity  deviations  in  the  9-state 
estimator,  only  the  7-state  estimator  is  considered  for 
further  testing.  In  other  words,  the  jerk  states  are  no 
longer  modelled,  and  only  the  7-state  propagation  equations. 
Equations  (34)  to  (39),  hold  henceforth. 

To  recap  the  findings  of  the  above  test,  it  is  concluded 
that,  with  the  current  values  of  R and  2.  dropping  the  jerk 
states  produces  better  estimator  performance.  A white  noise 
is  an  adequate  model  of  jerk,  but  Brownian  motion  is  not.  A 
suitable  representation  of  the  aircraft  dynamics  is 

White 

Noise 
> 

Suggestions  for  investigating  other  models  are  also  reserved 
f pr  the  recommendations . 

Comparison  of  Estimator  and  FASTtlAP  Filter  Positions  and  Velocities 

Figure  9 compares  FASTMAP  and  estimator  trajectories  for 
the  first  33  seconds  of  the  test  flight.  Due  to  the  identical 
'initial  conditions  both  begin  in  nearly  perfect  agreement  for 
the  first  three  seconds.  After  the  estimator  starts  depending 
on  its  own  propagation  and  update  equations,  the  estimator 
and  the  FASTMAP  filter  diverge  for  this  time  period  to  a 
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maximum  separation  of  180  feet  at  t equal  to  14.5  seconds. 
(This  information  is  obtained  from  output  data) . 

An  important  note  to  emphasize  is  that  the  FASTMAP 
filter  position  does  not  represent  true  position.  In  other 
words,  the  difference  vector  d^p 


-EF  -estimator  " -FASTMAP 


is  not  meant  to  be  interpreted  directly  as  an  error  e^,  : 


e = X - X 

— E —estimator  —true 


Rather,  the  estimator  performance  is  just  compared  to  the 
FASTMAP  filter.  Table  V illustrates  the  differences  in  feet 
between  the  filter  and  the  estimator  indicated  position.  The 
differences  in  the  last  column  are  computed  using  the 
Pythagorean  theorem: 

Difference  = (Let.  Difference')^  ^ CLong.  Difference) 
\ in  feet  / \ in  feet  / 

(73) 

Table  V.  Comparison  of  FASTMAP  Filter  and 
; Estimator  Position 


Time 

(sec) 

No . of  Measurements 
Since  Initiation  of 
Test  Flight 

Distance  Between 
Estimator  and 

Filter 

3.43 

20 

33.25 

6.12 

40 

135.34 

10.22 

60 

80.03 

14.41 

80 

178.00 

18.91 

100 

80.08 

22.97 

120 

57.60 

26.97 

140 

61.31 

31.34 

160 

74.99 
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Figure  11.  FASTMAP  and  Estimator  Longitude  for  First  33  Seconds 


Figures  10  and  11  show  the  filter  and  estimator  latitude 
versus  time,  and  FASTMAP  and  estimator  longitude  versus  time, 
respectively.  The  larger  deviations  occur  in  longitude, 
while  latitude  differences  between  the  estimator  and  FASTMAP 
remain  small. 

An  explanation  for  the  larger  longitude  differences  is 
obtained  from  a study  of  the  geometry  of  the  stations  relative 
to  the  initial  aircraft  position.  If  all  the  stations  lie  in 
the  same  general  direction  relative  to  the  aircraft,  then 
accuracy  in  that  direction  can  well  be  expected  to  be  better 
than  in  the  other  orthogonal  direction.  During  the  time  that 
Figures  10-11  cover,  stations  with  channel  numbers  57,  37,  59, 
and  43  are  acquired.  The  geometry  of  these  four  stations 
relative  to  the  aircraft  is  shown  in  Figure  12. 
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Figure  12.  Relative  Station  Geometry  For 
First  33  Seconds  of  Flight 


From  Figure  12,  it  is  noted  that  the  four  stations  supply 
ranges  in  a more  East  - West  direction  than  in  a North  - South 
direction.  This  effectively  causes  the  estimator  to  achieve 
more  benefit  from  the  measurement  information  in  the  longitude 


direction  than  in  the  latitude  direction.  Since  better 


measurement  information  is  available  in  the  East  - West 
direction,  the  H matrix  has  a value  of  the  partial  of 
longitude  with  respect  to  range  greater  than  that  of  the 
latitude  partial  which  in  turn  produces  a larger  Kalman  gain 
for  the  longitude  (East  - West)  direction.  Thus,  the 
latitude  states  depend  on  the  system  dynamic  equations  more 
heavily  than  the  longitude  states.  Because  of  the  constant 
velocity  (zero  acceleration)  model  that  is  initially  assumed 
and  the  lower  latitude  Kalman  gains,  the  latitude  state 
estimates  show  a constant  velocity  trend  as  depicted  in  Figure 
10.  Since  FASTMAP  also  indicates  an  almost  constant  velocity 
trajectory  for  the  first  30  seconds  of  flight,  FASTMAP  and 
the  estimator  closely  agree  in  lattitude  position  estimates. 

Estimator  longitude  states  on  the  other  hand  rely  more  on 
the  measurement  history  rather  than  on  the  constant  velocity 
model  because  of  the  higher  Kalman  gains. 

The  main  point  to  be  drawn  from  Figures  10  and  11  and 
the  preceding  note  is  that  the  accuracy  of  the  position 

j 

estimates  in  both  directions  depends  on  the  geometry  of  the 
stations  with  respect  to  the  aircraft.  As  the  flight 
proceeds,  different  stations  are  acquired,  and,  therefore,  i 

the  geometry  constantly  changes.  (Also,  the  forward  motion 
of  the  aircraft  causes  the  geometry  to  change  continuously.) 

Changing  geometry  in  turn  causes  repeated  transitions  from  a j 


model  dependence  of  the  states  to  a measurement  dependence 
and  back,  depending  on  the  Kalman  gain  value  (K,,  ,,«.). 


FASTMAP  and  the  estimator  latitude  velocities  are 
compared  in  Figure  13.  Approximately  one  second  after  the 
test  flight  initialization,  the  velocities  begin  to  differ. 
However,  these  differences  are  not  biased  in  any  one  direction, 
but  rather  switch  signs  throughout  the  30  second  interval. 

Due  to  the  gap  between  FASTMAP  and  the  estimator 

longitude  position  in  Figure  11,  longitude  velocities  are 

anticipated  to  be  quite  different.  As  Figure  14  indicates, 

longitude  velocities  are  indeed  different.  In  fact,  the 

widest  margin  has  a velocity  difference  of  almost  25  nm/hr . 

(See  Figure  14.)  This  may  seem  to  be  wholly  unacceptable, 

but  three  points  are  to  be  noted.  First  the  estimator  is 

2 

not  tuned;  in  other  words,  R is  still  1000  ft  and  the  ^ 

-17  2 5 

elements  are  1 x 10  rad  /sec  . Second,  FASTMAP  is  only  a 
reference  and  not  an  absolute  truth  model.  Finally,  the 
figure  shows  that  the  two  filters  begin  at  different  initial 
conditions . 

The  justification  for  beginning  the  estimation  process 
at  a slightly  different  initial  velocity  is  to  illustrate 
the  corrective  nature  of  the  estimator.  Initial  FASTMAP 
longitude  velocity  is  95.73  nm/hr  in  a westerly  direction, 
and  initial  estimator  longitude  velocity  is  set  at  93 
nautical  miles  per  hour  in  the  same  direction.  In  the  first 
'12  to  13  seconds  the  two  velocities  diverge.  After  that 
time,  the  velocity  difference  decreases  until  t is  equal  to  19 
seconds,  where  the  difference  remains  near  zero.  The 
estimator  essentially  is  shown  by  Figure  14  to  be  able  to 
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Figure  13.  FASTMAP  and  Estimator  Latitude  Velocities  for  First  33  Seconds 
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Figure  14.  FASTMAP  and  Estimator  Longitude  Velocities  for  First  33  Seconds. 
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track  the  aircraft  velocity  to  some  degree  of  accuracy 
even  when  exact  initial  conditions  are  not  known. 


Estimator  State  Variances 


Latitude  position  variance  P^^^(l,  1),  longitude 
position  variance  ^) . and  bias  variance  7), 

are  plotted  for  the  first  record  of  data  in  Figure  15.  The 
purpose  of  this  chart  is  to  observe  the  relative  characteristics 
of  these  variances  during  the  first  record,  and  not  to  study 
the  actual  variance  magnitudes . The  values  on  the  vertical 
axis  represent  the  variance  magnitudes  with  the  current  values 
for  R and  5,*  Since  R and  are  subject  to  change,  variance 
magnitudes  are  also  subject  to  change. 

Only  updated  variances,  P^j^^(l,  1),  and 

P(k)(7,  7)  are  plotted  in  Figure  15;  in  other  words,  the 

variance  is  plotted  only  after  each  measurement  is  taken.  The 

propagated  variances  are  not  shown  in  the  figure  so  as  to 

facilitate  the  viewing  of  the  overall  trend  in  the  variances . 

It  is  apparent  from  Figure  15  that  the  bias  variance  begins 

at  .75  X 10  rad^,  or  .25  x 10  rad^  less  than  the  initial 

-12  2 

variance  of  1 x 10  rad  . This  is  a direct  consequence  of 
plotting  only  the  updated  variances . The  variances  at  time 
equal  to  zero  are  plotted  after  the  first  measurement  is 
taken,  which  are,  as  expected,  lower  than  the  initial 
variances . 

Generally,  the  bias  variance  curve  decreases  over  the 
6.29  seconds  interval.  However,  the  appearance  of  two  sharp 
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peaks  in  the  curve  indicates  sudden  increases  in  bias  uncer- 
tainty. The  first  peak  occurs  at  t = .99  seconds  when  station 
43  is  called  on  for  the  first  time,  and  the  second  peak 
occurs  when  channel  59  is  acquired  for  the  first  time  at 
t = 3.12  seconds.  When  a station  is  called  on  for  the  first 

time,  bias  for  that  station  is  completely  unknown;  bias 

2 

uncertainty  is  reinitialized  to  a = 432  ft  at  each  new 
station  acquisition,  causing  the  peaks  in  the  bias  variance 
curve  of  Figure  15 . 

Latitude  and  longitude  position  variances  are  represented 
by  the  other  two  curves  of  Figure  15 . They  begin  at  the  low 
initial  value  and  grow  throughout  the  first  record.  To  show 
that  the  position  variances  actually  do  reach  an  upper  bound. 
Figure  16  carries  latitude  variance  through  the  first  four 

^ 4- 

records  of  data.  At  t ==  10  seconds,  P^j^^(l,  1)  halts  its 
ascent  and  settles  to  a value  somewhat  lower  than  the  peak. 

2 

At  t = 166.9  seconds  latitude,  variance  is  equal  to  105.91  ft 
or  2.45x10”^^  rad^ . (See  Figure  16.) 

The  fact  that  station  bias  values  are  estimated  with 
better  accuracy  as  more  measurements  arrive  can  be  a possible 
explanation  for  the  variance  overshoot  in  Figure  16.  DME 
measurements  are  predicted  with  increasing  accuracy  as  time 
goes  from  t = 10  seconds  to  the  next  station  acquisition 
due  to  ever-increasing  bias  knowledge.  With  decreasing 
measurement  uncertainty,  total  system  variance  (or  uncertainty) 
also  drops,  creating  the  gradual  decrease  in  the  variance 
plot . 


64 


10.00 


DME  Station  Bias  Results 


Figxire  17  shows  the  estimator  bias  values  for  the  four 
stations  used  in  the  first  record.  Channels  57  and  37  come 
in  at  t = 0 seconds  and  t = .17  seconds  respectively,  channel 
43  comes  in  at  t = .99  seconds,  and  channel  59  comes  in  at 
t = 3.12  seconds.  I'Jhen  each  station  is  called  on  for  the 
first  time,  the  bias  prediction  changes  the  most,  but  the 
more  measurements  that  are  used,  the  less  the  bias  values 
change.  This  is  a direct  consequence  of  the  gradual  decreasing 
bias  variance  curve  of  Figure  15.  In  other  words,  bias 
predictions  must  change  less  for  lower  uncertainty. 

The  most  obvious  characteristic  of  the  bias  curves 
in  Figure  17  is  the  tendency  to  flatten  out  as  more  measure- 
ments are  used.  The  final  bias  values  for  the  first  four 
stations  during  the  first  record  can  be  summarized  as 

Channel  57  =>  6.71  feet 

Channel  37  =>  -313.16  feet 

Channel  43  =>  -74.13  feet 

Channel  59  =>  -257.60  feet 

Of  course,  these  values  are  also  subject  to  change  as  the 
flight  progresses  because  of  changing  geometry. 

To  illustrate  that  the  residuals  in  general  actually  do 
decrease  as  better  bias  estimates  are  available  for  the 
prediction  equation.  Figures  18  and  19  show  the  residual  and 
bias  history  of  stations  37  and  59  for  the  first  record  of 
data.  Both  station  biases  approach  an  apparent  steady  state 
value,  while  the  residuals  (those  corresponding  to  the  same 
stations  and  the  same  time)  generally  decrease. 
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Figure  17.  Four  Station  Bias  Values  for  First  Record  Showing  A Steady  State  Tendency 
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During  the  course  of  the  test  flight,  many  additional 

stations  are  acquired.  Approximate  bias  steady  state  values  j' 

1 

[ are  obtained  from  the  estimator  for  all  the  DME  stations  ! 

\ 4 '' 

I during  the  first  40  minutes  of  the  test  flight.  Since  the  | 

I FASTMAP  filter  does  not  calculate  station  bias  estimates,  the  ! 

; ! 

accuracy  of  the  estimator  bias  estimates  is  not  evaluated  ' 

\ 

in  this  study.  Therefore,  these  estimates  are  only  presented 
in  Table  VI  to  show  the  results  of  the  bias  model.  (See 
Chapter  III.)  Each  approximate  steady-state  bias  quantity 
in  the  table  represents  the  last  estimated  value  on  or  before 
the  40  minute  mark . 


Table  VI.  Summary  of  Approximate 
DME  Stations  Used  During  First 

Steady  State  Bias  for 
^0  Minutes 

STATION  NAME 

IDENTI- 

FICATION 

CHANNEL 

APPROXIMATE 
STEADY  STATE 
BIAS  ESTIMA- 
TION (FEET) 

Dalhart 

DHT 

57 

14.5 

Anton  Chico 

ACH 

37 

-1.7 

Santa  Fe 

SAF 

43 

207.9 

Texico 

TXO 

59 

-352.0 

Newman 

EWM 

71 

-289.2 

Truth  or  Consequences 

TCS 

74 

-267.1 

Alburquerque 

ABQ 

79 

119.1 

Holloman 

HMN 

85 

80.3 

Corona 

CNX 

102 

164.5 

Farmington 

FMN 

100 

192.7 

Roswell 

ROW 

108 

-434.5 

Cannon 

CVS 

104 

96.0 

Cimarron 

CIM 

111 

178.8 

Socorro 

ONM 

115 

-40.2 

Las  Vegas 

LVS 

120 

-65.0 

Taos 

TAS 

123 

-164.0 

Zuni 

ZUN 

81 

-331.7 

Tucumcari 

TCC 

83 

-277.7 

' Alamosa 

ALS 

86 

-310.4 

Colorado  Springs 

COS 

72 

-87.9 

Tobe 

TBE 

105 

-399.5 

Pueblo 

PUB 

114 

-71.8 

Thurman 

TXC 

76 

407.6 

Denver 

DEN 

110 

-91.6 

71 
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1 V.  Estimator  Performance  Compared  to  CIRIS 
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CIRIS  Pata 

In  Chapter  IV  it  is  shown  that  the  designed  estimator 
compares  well  to  the  FASTMAP  filter  (an  accepted  and  usable 
filter).  However,  since  FASTK/iP  states  cannot  be  considered 
as  true  states,  the  FASTMAP  filter  only  offers  an  approximate 
method  for  testing  the  cotimatcr.  A better  method  of  testing 
is  now  sought  to  provide  a more  acceptable  error  analysis  for 
the  estimator. 

Comparison  to  CIGTF's  Completely  Integrated  Reference 
Instrumentation  System  (CIRIS)  that  recorded  trajectory  data 
(along  with  the  FASTI-IAP  filter)  for  the  same  test  flight  is 
the  approach  used  in  this  chapter.  CIRIS  calculated  verti- 
cal, longitude,  and  latitude  position  and  velocity  for  the 
entirety  of  the  test  flight.  These  six  states,  along  with  a 
standard  time  (IRIG  time)  and  64  other  words  of  information, 
are  all  stored  for  each  set  of  calculations  (one  file)  on  a 
9- track  tape. 

CIRIS  Estimates  of  Error 

Although  a comparison  to  CIRIS  states  can  be  considered 
a more  accurate  method  of  testing  the  estimator  than  a com- 
parison to  the  FASTMAP  filter  states,  CIRIS  also  has  asso- 
ciated errors.  An  ideal  evaluation  of  estimator  performance 
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would  necessitate  the  use  of  the  error 


-E 


-ESTIMATOR 


X 

“TRUE 


(74) 


However,  a comparison  of  estimator  states  to  CIRIS  states 
yields  results  characterized  by 


-ESTIMATOR  " “CIRIS 


(75) 


Therefore,  CIRIS  errors  must  be  noted,  which  are  represented 
by 


e 

-c 


= X 


-CIRIS 


-TRUE 


(76) 


From  Eqs  (74),  (75),  and  (76),  the  desired  result  can  then 
be  obtained  as 


(77) 


Assuming  zero-mean  errors  and  since  estimator  and  CIRIS 
errors  are  independent,  the  total  error  correlation  can  be 
expressed  as 


TOTAL  ERROR  r T, 
CORRELATION  " 


l-EC-ECJ 


+ E 


[e  e^^ 
L-C-i 


C J 


(78) 


'CIRIS  estimates  of  the  square  root  of  the  diagonal  terms  of 
the  last  term  in  Eq  (78)  are  located  in  words  35  through  40 
of  each  file.  These  errors  are  illustrated  in  Table  VII, 


Table  VII 


\ 


CIRIS  Estimates  of  Standard  Deviation  of  Error 


Quantity 


CIRIS  Estimate  of  Standard 
Deviation  of  Error 


^Longitude 


Latitude 


Altitude 


VE 


VN 


W 


Longitude  Standard  Deviation 
Latitude  Standard  Deviation 
Altitude  Standard  Deviation 
East  Velocity  Standard  Deviation 
North  Velocity  Standard  Deviation 
Vertical  Velocity  Standard  Deviation 


Table  VIII  shows  some  typical  values  of  CIRIS  errors  sampled 
about  four  minutes  after  the  initialization  of  the  FASTMAP 
test.  Time  corresponds  to  the  time  following  the  FASTMAP 
initialization.  With  CIRIS  errors  in  mind,  information  from 
the  CIRIS  tape  can  be  used  to  evaluate  estimator  performance. 

Matching  CIRIS  and  Estimator  States  in  Time 

IRIG  time  stored  on  the  CIRIS  tape  is  synchronous  to  the 
IRIG  time  recorded  on  the  FASTMAP  filter  tape.  The  estimator 
uses  IRIG  time  from  the  FASTMAP  filter  tape,  and  therefore  a 
matchup  between  CIRIS  and  estimator  states  is  conceivable. 
However,  both  CIRIS  and  the  FASTMAP  filter  have  their  own 
discrete  values  of  IRIG  time,  complicating  the  match-up 


L.  _ : 
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Table  VIII 


r 


Example 

of  CIRIS  Estimates  of  Error 

Time 

(sec) 

°Long . 
(feet) 

'^Lat.  ®VE 

(feet)  (ft/sec) 

°VN 

(ft/ sec) 

240.78 

107.4 

23.2 

.6900 

.1407 

247.86 

55.9 

28.3 

.4584 

.3102 

256.14 

39.5 

19.7 

.4380 

.2310 

procedure.  In  order  to  compare  properly,  some  method  is 
necessary  to  match  the  appropriate  CIRIS  states  to  the 
corresponding  estimator  states. 

Because  IRIG  time  is  only  available  at  different  dis- 
crete values  for  each  tape,  estimator  states  are  modified 
in  order  to  be  comparable  to  the  CIRIS  states  that  are 
closest  in  time  to  the  estimator  states.  Figure  20  illus- 
trates the  fact  that  the  estimator  and  CIRIS  have  mismatch'.d 
state  values.  In  the  figure,  the  term  "desired  estimator 
St, ate  values"  refers  to  the  estimator  states  that  directly 
correspond  in  time  to  the  available  CIRIS  states. 

To  obtain  the  desired  estimator  state  values,  state 
values  are  propagated  and  updated  in  the  manner  described  in 
Chapter  III  up  to  and  including  measurements  (i).  Next,  the 
time  at  measurement  (i)  is  subtracted  from  the  time  corre- 
sponding to  the  available  CIRIS  states.  The  time  difference. 
At',  is  employed  in  the  estimator  propagation  equations  to 
propagate  the  states  from  measurement  (i)  to  the  desired 
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Figure  20.  CIRIS  and  Estimator  States  Mismatched  in  Time 


estimator  state  values.  The  two  systems  now  have  synchronous 
state  values  that  match  in  time  and  are  directly  comparable. 
The  sample  time  is  7 seconds  and  total  number  of  sample 
points  are  400.  In  reference  to  Figure  20,  this  procedure 
can  be  summarized  in  four  steps  as  follows: 

STEP  1 : Propagate  and  update  states  to  measurement  (i) . 

STEP  2:  Obtain  the  time  difference  = At'  = CIPx.IS - 

TIME 

FASTMAP,j, , where  corresponds  to 

the  DME  measurement  (i)  just  before  the 
available  CIRIS  state. 

STEP  3 : Use  At'  in  the  estimator  propagation  equations 
to  propagate  the  states  to  the  desired  time. 

, STEP  4 : Compare  CIRIS  and  estimator  states. 

Test  Results 

Using  the  matching  procedure  of  the  previous  section, 
this  test  entails  the  comparison  of  the  estimator  and  CIRIS 


latitude,  longitude,  latitude  velocity  longitude 

velocity  . The  study  includes  the  test  flight  from 

the  start  of  incoming  DME  measurements  (the  investigation  of 
the  FASTMAP  filter)  through  the  180°  turn  and  about  8 minutes 
past  the  turn,  totaling  50  minutes  of  flight.  The  estimator 
test  is  extended  beyond  the  turn  to  show  that  the  estimator 
follows  the  turn. 

The  180°  turn  provides  a rigorous  test  for  the  estimator 
by  checking  its  performance  in  such  trajectory  chano-es . 

Slight  trajectory  changes  also  exist  in  the  test  flight  prior 
to  the  180°  turn.  Figure  21  shows  CIRIS  latitude  versus  time 
for  18  minutes  following  FASTMAP  initialization.  A slight 
increase  of  latitude  velocity  (in  other  words,  a small 
increase  in  slope)  is  indicated  at  approximately  t = 560 
seconds.  A more  prominent  change  in  trajectory  is  noted  in 
Figure  22,  a representation  of  CIRIS  longitude  versus  time. 
Longitude  velocity  changes  from  a westward  direction  (indi- 
cated by  the  decreasing  slope)  to  an  eastward  direction  at 
t ■=  560  seconds.  Velocity  estimates  from  the  estimator 
should  indicate  these  changes  in  trajectory. 

The  trajectory  changes  indicated  above  result  of  course 
from  changes  in  velocity.  Acceptable  estimator  performance 
.means  that  estimator  velocities  track  these  CIRIS  velocities. 
Figure  23  compares  both  latitude  and  longitude  velocities 
for  CIRIS  and  the  estimator.  The  smoother  upper  and  lower 
curves  represent  CIRIS  latitude  and  longitude  velocities 
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Figure  21.  CIRIS  Latitude  versus  Time  for  18  Minutes  of  Flight 


Figure  22.  CIRIS  Longitude  versus  Time  for  18  Minutes  of  Flight 


respectively,  and  the  curves  that  oscillate  about  CIRIS 
velocities  are  the  estimator's  own  latitude  and  longitude 
velocities.  As  Figure  23  shows,  the  estimator  velocity 
estimates  actually  do  track  the  changes  in  CIRIS  velocity. 

As  noted  in  Chapter  I,  high  frequency  errors  in  the  reference 
system,  such  as  the  estimator  velocity  curves  in  Figure  23, 
are  not  as  significant  as  unbounded  error  growth.  On  this 
basis,  small  estimator  state  oscillations  about  the  smoother 
CIRIS  states  are  acceptable.  (Note  that  tuning  could  remove 
them  to  some  degree.)  The  error  using  Eq  (75)  is  plotted 
versus  time  for  the  first  50  minutes  of  flight  in  Figures  24 
through  27  for  positions  and  velocities.  The  sharp  peaks  in 
all  the  error  curves  represent  the  180°  turn.  As  can  be 
expected,  estimator  accuracy  decreases  momentarily  for  the 
duration  of  the  turn.  However,  these  figures  show  that  the 
estimates  quickly  recuperate  from  the  turn.  Although  not 
attempted  in  this  study,  ^ can  be  made  to  increase  when  a 
turn  is  indicated  to  put  more  estimator  dependence  on  DME 
measurements  and  less  on  the  internal  model.  (See  recommen- 
dations . ) 

Figures  28  through  31  are  histograms  that  represent  the 
frequency  distribution  of  the  latitude  and  longitude  position 
and  velocity  errors.  Generation  of  these  histograms  is 
accomplished  by  dividing  the  error  (Eq  (75))  into  50  feet 
sections  for  position  and  5 feet  per  second  sections  for 
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Figure  28.  Latitude  Position  Error  Histogram 
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Figure  30.  Latitude  Velocity  Error  Histogram 
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Figure  31.  Longitude  Velocity  Error  Histogram 
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velocity,  and  noting  the  number  of  errors  between  each 
section;  error  'sampling  is  one  sample  every  seven  seconds. 

Figures  25  and  29  show  that  longitude  position  error 
is  biased  about  150  to  200  feet.  The  fact  that  bias  errors 
are  commonly  found  in  EKF  designs  may  be  a possible  expla- 
nation. Bias  errors  in  EKF  designs  can  result  from 
assumptions  of  the  form 

Z f(x)  (79) 

for  nonlinear  f,  such  as  in  the  calculation  of  the  linearized 
measurement  equation  in  Chapter  III.  Partial  corrections  for 
biases  in  EKF  designs  can  be  accomplished  by  using  "Bias 
Correction  Terms."  Use  of  these  terms  partially  corrects 
for  the  bias  by  adding  higher  order  terms  to  both  the  state 
propagation  and  the  measurement  prediction  equations.  This 
procedure  is  reserved  for  a possible  future  study  and  is  not 
attempted  in  this  report. 

Using  the  histograms  in  Figures  28  through  31,  approxi- 
mate values  can  be  obtained  that  encompass  50  percent  of  all 
the  errors  for  latitude,  longitude,  latitude  velocity,  and 
longitude  velocity.  The  magnitudes  of  these  values  should 
give  a rough  estimate  of  estimator  performance  and  feasibil- 
ity for  using  the  estimator  as  an  INS  reference  system.  The 
approximate  values  that  contain  50  percent  of  the  errors  are 
obtained  as 
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Latitude  position  + 83  feet 

Longitude  position  =>.  + 183  feet 

Latitude  velocity  + 8.4  feet/sec 

Longitude  velocity  =>  + 7.5  feet/sec 

These  approximate  values  exceed  the  accuracy  goals  stated 
in  Chapter  I.  However,  with  proper  tuning  and  with  the 
"Bias  Correction  Terms,"  the  accuracy  goals  might  be 
achievable . 


VI . Effects  of  Omitting  Measurements 


Omission  of  Bad  Measurements 

In  designing  the  estimator,  some  method  of  omitting 
erroneous  measurements  should  be  considered.  Estimator 
performance  depends  on  the  measurement  history.  If  one 
measurement  in  a group  of  measurements  is  obviously  different, 
then  a possible  action  would  be  to  rid  the  system  of  that 
measurement  since  keeping  it  might  well  degrade  performance. 

On  the  other  hand,  removing  a measurement  is  like  throwing 
away  a piece  of  informat ion- -information  which,  no  mater  how 
erroneous,  might  actually  help  the  estimation  process.  How 
bad  must  a measurement  be  before  no  helpful  information  can 
be  derived  from  it? 

One  way  to  define  "erroneous"  measurements  is  to  include 

all  measurements  with  residuals  beyond  an  na  „ boundary, 

Rho 

where  the  residual  standard  deviation,  a _,  is  given  by 
, Rho 


°RES  [^(k+1)  ^(k+1)  -(k+1) 


and  n is  a value  on  the  order  of  2 or  3. 

Estimator  values  of  are  plotted  in  Figure  32  for  the 

first  nine  records  (about  60  seconds)  of  measurement  data. 

Another  way  to  define  "erroneous"  measurements  in 
relation  to  this  design  would  be  to  include  all  measurements 
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with  residuals  over  "X"  feet,  with  "X"  selected  a priori, 

A concise  answer  to  the  above  question  would  require  an 
appropriate  value  fir  X.  (This  value  could  then  be  compared 

to  the  time  history  of  in  Figure  32.)  In  order  to 

determine  a "ball-park"  value  for  X,  a study  presented  below 
contrasts  estimator  performance  with  four  values  of  X.  All 
measurements  with  residuals  greater  than  X feet  are  omitted. 
Estimator  performance  is  approximated  by  noting  the  time- 
averaged  difference  between  estimator  postion  and  the  FASTIiAP 
filter  position  for  each  value  of  X.  Table  IX  shows  four 
values  of  X and  the  corresponding  position  differences  for 
both  latitude  and  longitude. 


Table  IX 

Four  Values  of  X and  Corresponding  FASTMAP 
and  Estimator  Position  Differences 


X 

(feet) 

Latitude 

Difference 

(feet) 

Longitude 

Difference 

(feet) 

CO 

, (all  measurements 

accepted) 

66.9 

98.3 

600 

71.6 

100.8 

400 

255.9 

452.6 

300 

767.1 

1818.2 

i 

1 

i Table  IX  reveals  that  when  all  measurements  with  residuals 

I 

I above  300  feet  are  omitted  from  the  estimator  equations,  the 


time-averaged  position  differences  between  the  estimator  and 


il 


the  FASTMAP  filter  are  relatively  large.  A plot  of  estimator 
variance  for  latitude  position,  (1 , 1)  , is  shown  in 

Figure  33.  From  the  figure,  the  latitude  position  variance 
can  be  seen  to  diverge  when  X is  set  at  300  feet.  As  the 
value  of  X increases  to  400  feet,  the  variance  drops  sharply. 
At  X equal  to  600  feet,  the  variance  levels  off  close  to  the 
value  for  X equal  to  infinity  (when  all  measurements  are 
accepted) . The  difference  between  X equal  to  600  and  infin- 
ity is  almost  negligible  in  terms  of  estimator  performance 
(See  Table  IX) . 

The  results  of  this  study  reveal  that  measurements  with 
residuals  over  600  feet  can  effectively  be  dropped  without 
losing  accuracy.  In  reference  to  Table  IX  and  Figure  33, 
the  minimum  value  for  X is  600  feet.  A value  for  X to  be 
used  in  the  actual  estimator  is  chosen  at  X equal  to 
2500  feet.  This  ensures  that  no  bad  measurements,  such  as 
wrong  station  coordinates  and  gross  receiver- transmitter 
delays,  are  used  in  the  estimator  equations.  On  the  other 
hand,  choosing  X well  over  the  minimum  (X  = 600  feet)  will 
ensure  response  to  trajectory  changes  such  as  180°  turns. 
Again,  2 values  can  be  made  to  increase  to  meet  trajectory 
changes  (See  Recommendations). 

Now  that  all  measurements  with  residuals  over  2500  feet 
are  omitted,  a method  for  implementing  this  is  necessary. 

The  omission  of  bad  measurement  data  is  shown  in  the  esti- 
mator flow  chart.  Figure  7.  A measurement  prediction  is 
made  from  both  prior  measurement  and  state  information  and 
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is  immediately  followed  by  the  actual  measurement.  If  the 
difference  between  the  prediction  and  actual  measurement  is 
greater  than  2500  feet,  the  estimator  treats  that  measurement 
as  if  it  never  existed.  The  estimator  carries  the  propa- 
gation of  the  states  across  this  measurement  time  to  the 
next  measurement.  However,  a sequence  of  residuals  con- 
sistently over  2500  feet  should  not  be  omitted,  since  this 
phenomenon  could  well  represent  some  trajectory  change. 
Another  future  study  not  attempted  in  this  report  would  be 
to  monitor  the  residual  sizes  and  increase  ^ when  sequences 
of  res ‘.duals  are  over  a certain  amount  (i.e,,  an  adaptive 
filter) . 

Effects  of  Increasing  the  Time  Between  Measurements 

The  time  difference  between  each  station  acquisition 
for  this  test  flight  is  on  the  order  of  .05  to  ,5  seconds. 
Current  DME  digital  equipment  allows  station  acquisition 
at  such  high  rates.  However,  if  station  acquisition  is  only 
possible  at  a slower  rate,  estimator  performance  is  expected 
to  degrade.  An  additional  study  presented  in  this  chapter 
illustrates  system  performance  as  the  time  between  each 
measurement  is  lengthened.  Performance  may  not  be  degraded 
to  an  undesirable  degree  when  measurements  are  acquired  at  a 
slower  rate.  Perhaps  there  is  a point  where  any  additional 
measurements  in  a given  time  interval  will  not  add  a sub- 
stantial amount  of  information  to  the  estimation  process. 

If  there  is  such  a point,  then  the  plot  of  estimator 
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Figure  34.  Relationship  between  the  Number  of 
. Measurements  Omitted  and  Time-Averaged 
I Estimator  Latitude  Variance 
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VII.  Conclusions  and  Reconmiendations 

Conclusions 

Based  on  the  design  and  testing  of  the  estimator,  the 
following  conclusions  are  drawn: 

1.  An  extended  Kalman  filter  using  multiple  DME  range 
measurements  is  a feasible  reference  system  for  examining 
INS  low  frequency  errors.  Estimator  errors  in  the  high  fre- 
quency range  are  greater  than  those  of  an  INS,  but  errors  in 
the  DME-based  estimator  are  consistent  in  strength  and  do 
not  exhibit  an  unbounded  growth  as  typical  of  INS  errors. 

2.  With  the  ^ R used  in  the  estimator  tests,  the 
7-state  estimator  without  the  jerk  states  provides  better 
performance  than  the  9-state  estimator  with  the  jerk  states 
included. 

3.  For  the  estimator  in  this  study,  the  approximate 
values  that  encompass  50  percent  of  all  the  errors  for  lati- 
tude, longitude,  latitude  velocity,  and  longitude  velocity, 
as  compared  to  CIRIS,  are  summarized  below: 

Latitude  position  + 83  feet 

Longitude  position  + 183  feet 

' Latitude  velocity  + 8.4  feet/sec 

Longitude  velocity  — S + 7.5  feet/sec 

Although  not  attempted  in  this  study,  properly  tuned  values 
of  R and  2 might  result  in  lower  estimator  position  and 
velocity  errors. 
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4.  The  estimator  is  able  to  track  changing  trajectories 
such  as  a 180°  turn.  Errors  in  the  estimator  states  do 
momentarily  increase  for  such  a change,  but  the  estimator 
rapidly  recovers  following  the  trajectory  change. 

5.  Measurement  residuals  over  600  feet  (or  more)  can 
be  dropped  without  significant  loss  of  estimator  accuracy. 

On  the  other  hand,  the  time  between  measurements  cannot  be 
lengthened  without  loss  of  estimator  accuracy. 

6.  Finally,  the  estimator  recovers  from  initial  con- 
ditions that  are  not  precisely  known. 

Recommendations 

The  following  recommendations  are  suggested  for  improving 
the  performance  of  the  multiple  DME  estimator  and  for  further 
studies  in  this  area: 

1.  Proper  R and  ^ values  for  the  estimator  should  be 
obtained  by  a covariance  and/or  a Monte  Carlo  analysis  of 
performance. 

2.  An  estimator  using  an  adaptive  approach  should  be 
investigated;  that  is,  R and  Q can  be  made  functions  of 
time,  environment  parameters,  or  current  residual  values. 

For  example,  the  constant  R used  in  the  estimator  tests  can 
be  designed  to  change  in  relation  to  the  magnitude  of  the 
t)ME  range  measurements.  Also,  Q can  be  boosted  at  trajectory 
changes  either  by  monitoring  the  residuals  or  by  informing 
the  estimator  of  commanded  aircraft  trajectory  changes. 


3.  Although  the  7-state  estimator  provides  better  per- 


formance in  this  study,  other  approaches  to  designing  the 
estimator  are  suggested.  The  7-state  estimator  is  currently 
modeled  in  the  form 


White 


Noise  Acceleration  I Velocity  | Position 


T 

!‘ 

This  study  would  search  for  the  optimal  combination  of  the 
white  noise  strengths  for  both  noises  and  Parameters 

to  vary  in  the  preceding  diagram  would  include  R,  the 
strengths  of  the  white  noises,  and  T,  and  a. 


103 


m 


Bibliography 


1.  AFM  51-40.  Air  Navigation . Washington:  Department  of 
the  Air  Force  and  the  Navy,  July  1973. 

2.  Kay tor,  Myron  and  Walter  R.  Fried.  Avionics  Navigation 
Systems . New  York:  Wiley,  1969. 

3 . Krumm , J . A . , et  al . Flight  Testing  of  the  Griarnrnan 
FASTMAP  System.  Unpublished  report.  Holloman  Air  Force 
Base,  New  Mexico:  Central  Inertial  Guidance  Test 
Facility,  May  1977. 

4.  Flight  Testing  of  the  Sierra  SCORE  System. 

Unpublished  report.  Holloman  Air  Force  Base,  New  Mexico; 
Central  Inertial  Guidance  Test  Facility,  May  1977. 

5.  Latham,  R.  W.  "Aircraft  Positioning  with  Multiple  DME." 

Navigation:  Journal  of  the  Institute  of  Navigation,  21: 

150-158  (Summer  1974). 

6.  Latham,  R.  W.  and  R.  S.  Townes.  "DME  Errors,"  Navigation: 
Journal  of  the  Institute  of  Navigation,  22:  332-342 
(Winter  1975-76yT~^ 

7.  Maybeck,  P.  S.  Filter  Design  for  a Tacon-Aided  Boro- 
Inertial  System  with  ILS  Smoothing  Capability. 
AFFDL-TM-74-52 . Wright-Patterson  Air  Force  Base,  Ohio: 
Air  Force  Flight  Dynamics  Laboratory,  1974. 

8.  McCaskill,  T.,  et  al.  "A  Sequential  Range  Navigation 

> Algorithm  for  a Medium  Altitude  Navigation  Satellite," 
Navigation:  Journal  of  the  Institute  of  Navigation,  23 : 

166  (Summer  1976) . 


1 


104 


J 


Appendix  A 


Local  Geodetic  Frame  to  XYZ  Coordinate  Transformation  (Ref  8:  166) 


A point  on  or  near  the  earth's  surface  can  be  described 
by  three  orthogonal  components,  X,  Y,  and  Z,  in  terms  of  its 
latitude,  longitude,  and  height  above  sea -level.  The  equations 
are  as  follows: 


1 

a 

4. 

1 

X =1 

-(l-e^sinV 

Hgt  . 

J cos  P. 

♦g 

cos  Pj 

(81) 

- a 

+ Hgt 

] \ 

sin  P;^ 

(82) 

Y = 1 

- (1-e  sin  Prf,  ) 

8 

- a 

■ + Hgt  - 

2 

ae 

"1 

sin  P 

8 

Z = 

(1-e^sin^P . 

'^8 

T,  2 . 2, 
(1-e  sin  . 

8 

where 

X 

Y 

Z 


a 

e 


Hgt 


Cartesian  coordinates  of  aircraft 

equatorial  radius  of  earth 
eccentricity  of  the  earth 
aircraft  latitude 

aircraft  longitude 
aircraft  height 


Appendix  B 


Linearized  Measurement  Equation 


The  non-linear  measurement  equation 


is  linearized  to  obtain 


(k+1) 


3c  = I ]«!>*  = I ] «P 

‘J’g  g ^ ^ 


-(k+1) 


-(k+1) 


+ [3c/3b.  1_  ]6b. 


-(k+1) 


and  thus  an  H matrix  becomes 


H I . 0,  0,  3e/8p  I 


g 


-(k+1) 


0,  0,  0,  0,  8e  3b, 


-(k+1) 


-(k+1) 


where 


= aircraft  latitude 


= aircraft  longitude 


*i(k+l)  = measurement  bias  prediction 

X =1  Cartesian  coordinates  of  aircraft 
u 

„ I using  formulas  in  Appendix  A 
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Cartesian  coordinates  of  the 
DME  stations  using  the  formulas 
in  Appendix  A 


= aircraft  altitude 


Expressing  equation  (84)  in  terms  of  P,  , , and  Hgt  yields 


■[(([ 
. (ft 


(1-e  sin  P,  )■ 


(1-e^sin^P,  ) 


(l-e^sin  P,  ) 
<P 

g 


Hgt  ] 


cos  P,  cos  P, 

(p  A 

g 


Hgt  J cos  P,  sin  P^ 


1-3) 

'"s) 


(1-e  sin  P,  ) 
■*g 


tm  ■ 


(k+1) 


The  partial  of  Equation  (86)  with  respect  to  latitude  becomes 


X [2(X^-X^)  3X^/9, 


+ 2(Y  -Y  ) 3Y  /8„ 
'us  u P, 


'(k+1) 


+ 2(Z  -Z_)  9Z  ./9„  I 1 


First  3X  /3P,  is  determined; 

u 4,g 


3X  /3„ 
u P , 


(1-e'^sin'^P 


-^  + Hgt  ) 


cos  P.  cos  P, 


] (88: 


3X  /3P, 
u (|) 


= -( T-^ F-  + Hgt) 

V(l-e^sin^P,  ^ 


cos  P,  sin  P, 

X <p. 


+ ae^ (1-e^sin^P , ) sin  P,  cos^  P,  cos  P^^ 

‘’’g  ^g 


Similarly  3Y^/3P^  is  determined  as  follows-. 


f /3P  = -f 

u'  g \ 


(l-e^sin'^P . 


+ Hgt)  si 


in  P,  sin  P. 

X <]), 


+ ae^d-e^sin^P.  sin  P,  cos^  P.  sin 

g g g (90 


and  3Z  /3Px  is 
g 


3Z  /3P  = (l-e^)  ae^  (1-e^sinV  sin^  P.  cos  P. 

^g  g g g 

2 

+ r 7 — T y-  2~^ ¥ J ^ 

^•■(l-e  sin^P.,  (1-e^sin^P^ 


(1-e  sin  P. 


Substituting  Equations  (89),  (90),  and  (91)  into 


Equation  (87)  yields  an  expression  for  8e/8P.  . The  values 

for  Pa  and  P,  are  obtained  from  the  estimator's  best  estimates 
<pg  A 

of  the  states.  3e/  9Pi  is  obtained  in  a similar  manner 


3e/9P,  = r%  (X  -X  )^  + (Y  -Y  + (Z  -Z  )^  1"^ 

X |_^'us  us  ^us  j 


x[2(Xu-Xp 


+ 2(Y  -Y  ) 3Y  /3P- 
'us'  u X 


-(k+1) 


+ 2(Z^-Zp  8Z^/DP^ 


-(k+1) 


■(k+1) 


9X^/3 P^ 


(1-e  sin  P, 


+ Hgt]  cos  P^  ) sin  P^  (93) 


9Y  /9P, 
u X 


(1-e  sin  P. 


+ Hgt]  cos  P^  ) cos  P^  (94) 


9Z  /3^  = 0 

u'  P^ 


Substituting  Equations  (93),  (94),  and  (95)  into  Equation 
(92)  yields  an  expression  for  9£/3P^.  The  perturbation  of  the 
nonlinear  measurement  equation  (Equation  (84))  is  now  in  the 
form  of  a linear  function  - an  H matrix.  The  units  of  the 
H matrix  components  are  ft /radian. 
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"^This  report  is  directed  toward  the  design  of  a real-time 
estimation  algorithm,  a Kalman  filter,  that  estimates  aircraft 
position  and  velocity  using  multiple  DiVS  range  measurements.  The 
estimator  is  designed  and  tested  for  feasibility  as  a I'eference 
system  for  examining  Inertial  Navigational  System  (INS)  lo-/ 
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Vifith  the  tuninfF  parameters  used  in  the  estimator  tests,  the 
?-Etate  estimator  provides  better  performance  than  the  9“State 
estimator.  An  approximate  analysis  of  the  7-state  estimator 
performance  (by  comparison  to  FASTMAP,  a currently  used  and 
accepted  filter,  and  CIRIS,  the  Completely  Integrated  Reference 
Instrumentation  System, ) reveals  that  estimator  errors  in  the  high 
frequency  range  are  greater  than  those  of  an  INS,  but  errors  in 
the  DMii-based  estimator  are  consistent  in  strength  and  do  not 
exhibit  an  unbounded  growth  as  typical  of  INS  errox's.  For  the 
estimator  in  this  study,  the  approximate  values  that  encompass  50 
percent  of  all  the  errors  (as  compared  to  GIHI3)  for  latitude, 
longitude,  latitude  velocity,  and  longitude  velocity  wei'e 


Latitude  position 
Longitude  position 
Latitude  velocity 
Longitude  velocity 


+ 83  feet 

+ 183  feet 
+ 8,4  feet/second 

i 7,5  feet/second 


Improving  estimator  performance  is  suggested  by  proper  tuning 
and  by  using  an  adaptive  approach. 
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